Satellite Remote Sensing of Global Land Surface Temperature: Definition, Methods, Products, and Applications

Land surface temperature (LST) is a crucial parameter that reflects land–atmosphere interaction and has thus attracted wide interest from geoscientists. Owing to the rapid development of Earth observation technologies, remotely sensed LST is playing an increasingly essential role in various fields. This review aims to summarize the progress in LST estimation algorithms and accelerate its further applications. Thus, we briefly review the most‐used thermal infrared (TIR) LST estimation algorithms. More importantly, this review provides a comprehensive collection of the widely used TIR‐based LST products and offers important insights into the uncertainties in these products with respect to different land cover conditions via a systematic intercomparison analysis of several representative products. In addition to the discussion on product accuracy, we address problems related to the spatial discontinuity, spatiotemporal incomparability, and short time span of current LST products by introducing the most effective methods. With the aim of overcoming these challenges in available LST products, much progress has been made in developing spatiotemporal seamless LST data, which significantly promotes the successful applications of these products in the field of surface evapotranspiration and soil moisture estimation, agriculture drought monitoring, thermal environment monitoring, thermal anomaly monitoring, and climate change. Overall, this review encompasses the most recent advances in TIR‐based LST and the state‐of‐the‐art of applications of LST products at various spatial and temporal scales, identifies critical further research needs and directions to advance and optimize retrieval methods, and promotes the application of LST to improve the understanding of surface thermal dynamics and exchanges.


Introduction
Land surface temperature (LST) is an important variable within the Earth's climate system (Dash et al., 2002; Z.-L. Li et al., 2013). It is an indicator of the energy and water exchange between the land surface and atmosphere, influences the rate and timing of plant growth (Myneni et al., 1997;S.-S. Peng et al., 2014;X. Zhang et al., 2007), and has been applied in various fields to study various subjects, for example, the surface energy balance (Friedl, 2002;Tajfar et al., 2020;Vancutsem et al., 2010;Xu, Guo, et al., 2019;Xu, He, et al., 2019), the urban heat island effect (C. Alexander, 2020;Fu & Weng, 2016;Masoudi & Tan, 2019;Menon et al., 2010;J. Peng et al., 2018;Wu et al., 2019), surface soil moisture (SM) and evapotranspiration (ET) estimations (Colliander et al., 2017;Gallego-Elvira et al., 2019;Kalma et al., 2008;Karnieli et al., 2010;R. L. Tang & Li, 2017a, 2017bW. Zhao et al., 2021), and climate change (Bechtel, 2015;Bright et al., 2017;Cai et al., 2018;Keenan & Riley, 2018;Y. Shen et al., 2021;Tran et al., 2017). New insights into LST could improve the understanding of the land surface-atmosphere exchange processes at global and regional scales and provide a valuable metric of surface state for various applications. Therefore, LST has been recognized as a key Earth surface parameter by the National Aeronautics and Space Administration, the United States (US) Climate Change Science Program, and as one of the 10 essential climate variables in the land biosphere designated by the Global Climate Observing System (Hollmann et al., 2013).
LST retrieval from TIR data has been garnering wide interest, and reviews on this topic have been conducted. For example, based on the NOAA Advanced Very High Resolution Radiometer (AVHRR) observations, Prata et al. (1995) provided a phased summary of LST estimation methods by focusing on discussions of atmospheric effects, geometrical effects, and emissivity effects. To complement the review in Prata et al. (1995), Dash et al. (2002) provided a more detailed review of the methods and their use of different satellite sensors. Moreover, Z.-L. Li et al. (2013) provided a more systematic and comprehensive review than those aforementioned by clarifying the assumptions, advantages, limitations, and requirements of LST retrieval and validation methods and presenting topics for further research. Nearly 10 years have passed since the latest review, and since then, new LST retrieval algorithms, newly launched TIR instruments, and new releases of LST products have been proposed. Furthermore, there is a rapid increment in the number of remotely sensed-based LST publications according to the searching results for journal articles indexed by the Web of Science (WoS) using the keywords "Land surface temperature" and "Remote sensing," and the publication number in 2020 is close to three times of that in 2013 ( Figure 1). Thus, there is a new era because of the increase in the number of TIR satellites supporting the rapid advances in LST retrieval and applications. Additionally, long-term LST products at different spatial resolutions are available, promoting the use of LST data as an important climate record. Many applications use TIR LST products, and users may experience problems, for example, having insufficient knowledge of the availability, uncertainties, and critical issues (e.g., spatial discontinuity and spatiotemporal incomparability) of those LST products. These problems prevent the development of TIR remote sensing and the effective use of LST products.
To clarify the aforementioned problems, this paper summarizes the development of TIR-based LST remote sensing and systematically reviews LST retrieval methods, available LST products, and major applications. This

Definition of Satellite-Derived LST
LST is the thermodynamic temperature of a thin layer in the interface between soil, vegetation or other surface components and the atmosphere. It reflects how hot or cold is the Earth's surface would feel to the touch. In the context of remote sensing, this thin layer (i.e., land surface) is a continuous projected surface comprising all visible components within the sensor instantaneous field of view (IFOV; Figure 2), and its thickness is of the order of penetration depth. The penetration depth is the distance after which the intensity of the radiation decreases by 1/e (about 37%), and can be expressed as (Hecht, 2002): where is the penetration depth, 0 is the wavelength used for the measurement, and is the imaginary part of the refractive index which is between 0.01 and 1. Since remote sensing LST is derived from radiometric measurements, for example, radiance or brightness temperature (BT), by eliminating atmospheric effects and correcting for the emissivity effect, it is called radiometric temperature (Norman & Becker, 1995). For a given satellite sensor, the radiometric temperature can be formulated as follows (Chandrasekhar, 1960;Dash et al., 2002;Prata et al., 1995): where s is the radiometric temperature; is the viewing zenith angle (VZA); is the viewing azimuth angle; λ is the channel-effective wavelength; −1 is the inverse function of Planck's law; , ↑ , and ↓ are at-sensor observed radiance, upward atmospheric radiance, and downward atmospheric radiance, respectively; is the channel atmospheric transmittance; is the channel land surface emissivity (LSE).
This radiometric temperature has four properties: 1. The definition of radiometric temperature (Equation 2) is scale-invariant, namely, it is applicable to any spatial scale (Becker & Li, 1995). 2. According to Equation 1, the penetration depth is about 0.1-10 times the wavelength used for the measurement. In TIR domain (the wavelength around 10 μm), the penetration depth approximately ranges from 1 to 100 μm. The depth is shallow, and this radiometric temperature is also called skin temperature (Norman & Becker, 1995). In microwave domain (the wavelength around 1 cm), the penetration depth approximately ranges from 0.1 to 10 cm, and this radiometric temperature is also called subsurface temperature Galantowicz et al., 2011;Grody & Weng, 2008). 3. Radiometric temperature is directional and dependent on the viewing angle. Thus, it is also referred to as directional radiometric temperature (Norman & Becker, 1995). 4. Radiometric temperature defined by Equation 2 represents an ensemble temperature from all independent homogeneous and isothermal components viewed within an IFOV (Norman & Becker, 1995): where T si , ε i , and a i are surface temperature, emissivity, and projected area weight for the ith homogeneous and isothermal component, respectively; and ∆R multi and ∆ε multi represent radiance and emissivity increments induced by the multiple scattering effect among the components, respectively. Notably, ∆R multi and ∆ε multi can be negligible if the surface is flat. In this case, surface emitted radiance is the sum of component surface radiances according to their area weights. Furthermore, radiometric temperature is approximately a linear combination of component temperatures weighted by the area if both component temperatures and emissivities have a small difference.
1. LST retrieval is underdetermined. If the radiance is measured in N TIR channels, there will be always N equations with N unknown LSEs plus an unknown LST, even if the atmosphere effect is removed. 2. LST retrieval is unstable. The measured radiances in N TIR channels are highly correlated, which makes the solution of the equations unstable.
To solve this ill-posed problem, additional assumptions and constraints are necessary to increase the number of equations or decrease the number of unknowns and decorrelate the set of equations consisting of different channels.
Since the 1970s, various LST retrieval algorithms have been developed and used to generate LST products, which can be classified into four categories: the single-channel (SC) algorithm, the split-window (SW)/dual-window (DW) algorithm, the temperature and emissivity separation (TES) algorithm, and the physics-based day/night (D/N) algorithm. The SC and SW algorithms require a priori knowledge of LSE; the TES and D/N algorithms require performing atmospheric correction in advance.

SC Algorithm
The SC algorithm is a method for LST retrieval implemented by inverting the radiative transfer equation (RTE) with the radiance measured in a single TIR channel within the atmospheric window, provided that the LSE is known in advance. The atmospheric quantities (transmittance, upwelling radiance, and downwelling radiance) involved in RTE are estimated using high-quality atmospheric transmittance/radiance code and accurate atmospheric profiles over the study area at the satellite overpass (Jiménez-Muñoz et al., 2009). Because of the difficulties in obtaining simultaneous atmospheric profiles, attempts have been made to alternatively estimate atmospheric quantities by using some empirical function of both total atmospheric water vapor content (WVC) and near-surface air temperature or of either (Jiménez-Muñoz & Sobrino, 2003;Jiménez-Muñoz et al., 2009;. The advantage of the SC algorithm is its simplicity and applicability to all TIR sensors due to the minimal TIR channel requirements. However, the limitations of the algorithm include that it (a) requires a priori knowledge of LSE of each pixel, (b) requires accurate atmospheric profiles or total atmospheric WVC or/and near-surface air temperature to estimate atmospheric quantities, (c) uses empirical relationships to provide poorer accuracies than using atmospheric transmittance/radiance code, and (d) must consider the effect of elevation on column water vapor content, air temperature, and other atmospheric factors.

SW/DW Algorithm
The SW algorithm retrieves the LST by using the differential atmospheric absorption in two adjacent TIR channels centered at 11 and 12 μm to correct for atmospheric effects, provided that the LSEs in both SW channels are known a priori (Becker, 1987;Becker & Li, 1990b;Wan & Dozier, 1996). To solve the underdetermined problem, atmospheric quantities were approximately parameterized as a function of the total WVC and the equivalent atmospheric temperature (T atm ), by consequence, provided that WVC could be estimated by using other methods, there are two channels' measured radiances with two unknowns, that is, LST and T atm . Using the first-order or second-order Taylor series expansion of the RTE, and eliminating T atm , various SW algorithms have been proposed to retrieve LST by using a generalized linear or non-linear form as follows (Z.-L. Li et al., 2013): = 0 + 1 + 2( − ) + 3( − ) 2 where T i and T j are the BTs measured in two SW TIR channels; b 0 -b 2 and c 0 -c 3 are coefficients that depend primarily on the spectral response function, that is, g i and g j , of the two channels, the two channel emissivities ε i and ε j , the WVC, and the VZA, or = ( , , , , WVC, VZA) ( = 0, 1, 2 or 3) The accuracy of the SW algorithm is dependent on the choice of the coefficients b k or c k . Different parameterization schemes of the coefficients b k or c k lead to different SW algorithms (Z.-L. Li et al., 2013).

Intercomparison Method
The intercomparison method compared two or more LST products to evaluate their consistency and uncertainty from their disparities (Adeniran et al., 2022;Jiménez et al., 2012;Qian et al., 2013;Sajib & Wang, 2020;Trigo et al., 2008Trigo et al., , 2021. This method is often used when there are no measurements of in situ LST, and synchronous LSE/atmospheric profiles, to evaluate the relative accuracy of LST products; it also serves as a satisfactory supplement to T-based validation (Trigo et al., 2008; and provides useful information regarding the spatial patterns of LST product discrepancies, including artificial artifacts caused by algorithmic issues, such as poor-quality auxiliary data (Duan, Li, Cheng, et al., 2017;B. H. Tang et al., 2015). Unlike the T-based method, the intercomparison method can manage a wide range of land surfaces and be applied over a long period and on a global scale. However, using this method alone is insufficient for comprehensive validation because various LST products could be consistent but biased with respect to the truth. Additionally, this method requires accurate match-ups in spatial resolution, observation angle, and acquisition time among LST products (Duan, Li, Cheng, et al., 2017;Ye et al., 2021).

Bibliometric Analysis
To access and evaluate the history and development of LST validation methods, a literature-based bibliometric analysis was conducted. In the Web of Science database, a query search using the phrases "land surface temperature" AND ("evaluation" OR "assessment" OR "comparison" OR "validation") retrieved 485 scientific articles, of which 99 records were selected and classified by considering which were the most relevant to LST validation methods (the data was collected on 27 October 2021). Figure 4 overviews a synthesis article citation network on the aforementioned three LST validation methods.
The network represents the research status and development of LST validation methods. For example, more articles are related to using the T-based method than the R-based and intercomparison methods. Among these articles, the highly cited articles or those published early usually represent the landmark developments in the LST validation bibliographic landscape. For instance, according to a highly cited reference published in 2002 by Wan et al. (2002), the T-based method began to emerge and mature since the MODIS payloads onboard Terra launched in 1999, although sporadic articles on LST validations existed before Wan et al. (2002). Subsequently, Jacob et al. (2004) first compared MODIS LST with ASTER LST by using the intercomparison method. In subsequent studies, the MODIS LST product was recognized as the most widely used well-validated LST reference data in intercomparison (Duan, Li, Cheng, et al., 2017;Gao et al., 2012;Qian et al., 2013;Zarei et al., 2021). Wan and Li (2008) published the first paper for purposing the R-based method. Afterward, refinements were made to the R-based method, providing a more comprehensive perspective of LST products Wan, 2014).

LST Product Levels
Due to the maturity of remote sensing technology, various LST products are gradually being released Ma et al., 2020;Wan, 2019).
To distinguish LST products and use them more appropriately, this study processed a series of LST products at various levels ranging from Levels 2 to 4, where a Level 2 LST product is generated from the radiance/BT product in Level 1 and lower-level LST products are the basis of higher-level LST products. The generation of higher-level products will introduce more spatial and temporal processing steps (Hulley et al., 2019). Figure 5 describes the logical relations of LST products for the levels. If the product level is given, you can judge which processes it has completed in the production process.
• A Level 2 product provides LST retrieved from a single orbit of a sensor at the same resolution. During the process, it may either remain in the latitude and longitude orientation as those of Level 1 radiance/BT products or undergo spatial projection processing for the user's convenience. Some quality control parameters (e.g., uncertainties, confidence flags, or cloud probability) may be included at this level. • A Level 3 product is derived based on a Level 2 product either from a single orbit of a single sensor or from multiple orbits of single/multiple sensors. A reprojection is conducted to resample this level product temporally or spatially into a gridded map projection format. • A Level 4 product is constructed by reanalyzing lower-level LST products from multiple sensors and combining them into a spatial/temporal grid to form a gap-free product to ensure completeness and consistency in providing a comprehensive science data output for Earth observing systems.

Typical LST Products
According to the differences in spatial resolution, LST products can be classified into three categories: high resolution, medium resolution, and low resolution. High spatial resolution is considered within the order of hundreds of meters, and medium spatial resolution and low spatial resolution are approximately the order of kilometers and several kilometers, respectively.
The typical high spatial resolution LST products are, for example, ASTER and Landsat LST products. The ASTER LST product (AST_08) provides surface temperatures of land areas acquired either during the daytime or nighttime at 90 m spatial resolution and every 16 days. It is generated along with the LSE product (AST_05) using the TES algorithm. The coordinate system of AST_08 is the Universal Transverse Mercator (UTM). The Landsat LST products include Collection 1 (Provisional) and Collection 2, which are produced with the SC algorithm under the provision of LSE from ASTER Global Emissivity Database. They are reorganized into the UTM projection for global or the Albers projection for the United States based on the WGS84 sphere at 30 m spatial resolution with the temporal resolution of 16 days. The Landsat and ASTER LST products are provided at Level 2. The MODIS series LST product (MxD11&21) is a typical product at medium spatial resolution. The term "MxD" refers to MOD and MYD, which represent the MODIS sensor onboard the EOS/Terra and Aqua satellites, respectively. Regarding the retrieval algorithm difference, the MxD series LST product can be separated into MxD11 and MxD21 and are provided with two product levels (Levels 2 and 3). MxD11_L2 is the daily 5-min Level 2 LST product with 1 km spatial resolution at the nadir direction. MxD11A1 and MxD11A2 are the daily and 8-day Level 3 1 km LST products, respectively. These LST products are produced by the SW algorithm. Correspondingly, MxD11B1, MxD11B2, and MxD11B3 are the daily, 8-day, and monthly Level 3 LST products at 6 km spatial resolution generated using the D/N algorithm. MxD11C1, MxD11C2, and MxD11C3 are the daily, 8-day, and monthly Level 3 LST products produced by resampling the corresponding MxD11B products into a 0.05° grid. Comparatively, the MxD21 products, including MxD21_L2, MxD21A/N, MxD21A2, and MxD21C1/C2/ C3, are generated based on the TES algorithm. MxD21_L2 is a Level 2 LST product at 1 km spatial resolution at nadir. MxD21A1D/N and MxD21A2 are the daily and 8-day Level 3 1-km products, respectively. MxD21C1/C2/ C3 are the daily, 8-day, and monthly Level 3 LST products produced at the 0.05° grid.
One of the typical low spatial resolution LST products is the Satellite Application Facility on Land Surface Analysis (LSA-SAF) product for MSG/SEVIRI and the LST product for Geostationary Operational Environmental Satellite-R (GOES-R)/Advanced Baseline Imager (ABI). The LSA-SAF LST product is processed using the SW algorithm, with the consideration of LST dependence on viewing illumination geometries (Ermida, Trigo, DaCamara, & Pires, 2018). It has a spatial resolution of 3-5 km and a temporal resolution of 15 min. The GOES-R ABI LST product is also generated using the SW algorithm, with the temporal resolution of 1 hr for every ABI Full Disk of the Earth, the Continental United States (CONUS) region, and the Mesoscale regions at the spatial resolution of 2 km in CONUS and Mesoscale and 10 km in Full Disk.
As shown in Table 2, most LST products are provided at Levels 2 or 3. By contrast, Level 4 LST products are relatively rare. The GlobTemperature project has released Level 4 LST products under the Data User Element of the European Space Agency's 4th Earth Observation Envelope Program. The Copernicus Global Land Service obtained from a constellation of Geostationary Earth Orbit (GEO) satellites and reliably provides the global land LST datasets (hourly LST and 10-day LST Daily Cycle), covering most of the globe's land surface within the −60°/70° latitude range at a spatial resolution of 0.05°. The hourly LST consists of instantaneous hourly LST for the areas covered by the GEO satellites. Meanwhile, the 10-day LST Daily Cycle provides the LST daily cycle over each 10-day period, including minimum, maximum and median LST values over the 10-days for each hour.
The main added-value of these two products is the description of extreme, median and the diurnal cycle LST values, which are not provided by other Low-Earth Orbit (LEO) satellites.
Generally, the higher the spatial resolution, the lower the temporal resolution. Landsat, the ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS), and ASTER LST products are the global products with the highest spatial resolution, and those of GOES, MSG, FY4, MTSAT/Himawari are the regional products with the highest temporal resolution. Regarding product naming, the rules are not unified, resulting in confusion for users. The map projections and retrieval algorithms also differ. In addition, the accuracy of some LST products is not claimed by the producer. All these factors make it difficult for users to select appropriate LST products. Therefore, an improvement would be to unify the LST product naming rule, from which the LST level and spatial/temporal resolutions can be easily identified. If possible, the same LST reference data should be used to validate those LST products so that users can select corresponding products according to the accuracy requirements of the application.

Intercomparison of LST Products
As shown in Table 2, more than 30 types of LST products are publicly available. Before applying those LST products, an understanding of their accuracy is necessary. Although there are many validations studies for different products using the T-based, R-based, and intercomparison methods, the validations are mainly conducted over limited areas and land surfaces and the product accuracy varied a lot even when using the same product over different land surfaces and at different time. Consequently, the uncertainty of the validation results may confuse the users of these LST products. To clarify this issue and achieve consistent accuracy in the assessments among different LST products, directly validating all LST products with the same reference data set is a prior option. However, there are insufficient worldwide high-quality in situ measurements over various land surfaces, resulting   in the T-based and R-based methods being infeasible. Consequently, the intercomparison method is an alternative way to assess the consistency among LST products.
To avoid the uncertainties induced by the differences in viewing location, local time, and geometry, four criteria were established to guarantee the reliability of the intercomparison results: 1. Collocation in space: The MxD11A1 product and the other 12 LST products are spatially collocated by aggregation or nearest-neighbor matching. Because of the relatively high confidence level in MODIS cloud masking, only MODIS pixels were flagged as the highest quality (i.e., quality control [QC] = 0) were used to identify valid pixel pairs. For fine spatial resolution LST products, for example, AST_08, its LST is aggregated to the MxD11A1 pixel scale by averaging all pixels within an MxD11A1 pixel. For medium spatial resolution LST products, for example, VNP21A1, a pixel nearest to the MxD11A1 pixel is selected. For coarse spatial resolution LST products, for example, ABI-L2-LST, MxD11A1 LST is aggregated to the coarse spatial resolution scale by averaging all MxD11A1 pixels within a coarse spatial resolution pixel. In this aggregation process, the percent of MxD11A1 pixels with QC = 0 within a coarse spatial resolution pixel is required to be higher than 50% considering the tradeoff between estimation accuracy and the number of valid pixel pairs. 2. Temporal concurrence: The difference in viewing time between two intercompared LST products is limited to be less than 15 min. 3. Viewing geometry alignment: The difference in VZA between two intercompared LST products is limited to be less than 15°. 4. Statistically significance: For each land cover type, only the results of each LST product with the number of valid pixel pairs larger than 2000 are considered statistically significant.
Notably, the Landsat and ECOSTRESS LST products are excluded because of a small amount of pixel pairs under the aforementioned criteria. Therefore, the intercomparison was only performed for the remaining 10 products.    croplands (CRP), cropland/natural vegetation mosaics (CNV), permanent wetlands (PWL), permanent snow and ice (PSI), water bodies (WTB), urban and built-up lands (UBL), and barren (BRN). To facilitate the analysis of the intercomparison results, these 17 land cover types were grouped into five categories: (a) forests (ENF, EBF, DNF, DBF, and MXF), (b) shrublands and savannas (CSR, OSR, WDS, and SVN), (c) grasslands and croplands (GRS, CRP, and CNV), (d) wetlands, snow and ice, and water bodies (PWL, PSI, and WTB), and (e) urban and barren surface (UBL and BRN). The results of the daytime comparison, in Figure 6a, are as follows: 1. The performance of the MxD11B1 product is highly consistent with that of the MxD11A1 product during daytime, with an absolute bias of less than 0.5 K and an RMSD of less than 1.5 K over all surface types. The reason is attributed to the fact that the MxD11A1 LST product was used as initial value to solve equations in the D/N algorithm which was used to generate the MxD11B1 product. Consequently, we do not further describe the intercomparison results of the MxD11B1 product over each specific surface type. 2. Forests. For ENF, except for SL_2_LST and EDLST, the other LST products are overestimated, especially MLST. Except for EDLST, ABI-L2-LST, and MLST, the other LST products have an RMSD of less than 2 K. For EBF, negative biases are obtained for EDLST and AGRI-LST. Only the RMSD of VLSTO is less than 2 K. For DNF, except for EDLST, the other LST products are overestimated. Only the RMSD of SL_2_LST is larger than 2 K. For DBF, negative biases are obtained for EDLST and AGRI-LST; positive biases are obtained for the other LST products, especially ABI-L2-LST. The RMSD of AST_08 and VLSTO is less than 2 K. For MXF, EDLST are underestimated. Except for EDLST, SL_2_LST, AGRI-LST, and MLST, the LST products have an RMSD of less than 2 K. 3. Shrublands and savannas. For CSR, except for AST_08, the other LST products have positive biases. VLSTO have an RMSD of less than 2 K. For OSR, except for EDLST, the other LST products are overestimated. Only AST_08 has an RMSD of less than 2 K. For WDS, except for EDLST and AGRI-LST, the other LST products have positive biases. Only VLSTO has an RMSD of less than 2 K. For SVN, only AGRI-LST are underestimated. Only AST_08 has an RMSD of less than 2 K. 4. Grasslands and croplands. For GRS, except for EDLST and AGRI-LST, the other LST products are overestimated, especially MLST. The RMSD of all LST products is larger than 2 K. For CRP, negative biases are obtained for EDLST and AGRI-LST; positive biases are obtained for the other LST products. Only ABI-L2-LST have an RMSD of less than 2 K. For CNV, except for EDLST and AGRI-LST, the LST products have positive biases, especially MxD21A1. Except for EDLST, VLSTO, and AGRI-LST, the RMSDs of the other LST products are larger than 2 K. 5. Wetlands, snow and ice, and water bodies. For PWL, all LST products have positive biases. Only VNP21A1 has an RMSD of less than 2 K. For PSI, except for AST_08 and EDLST, the other LST products are overestimated. Only AST_08 has an RMSD of larger than 2 K. For WTB, all LST products are overestimated, especially MLST. The smallest absolute biases are obtained for VLSTO. Except for MxD21A1, AGRI-LST, and MLST, the other LST products have an RMSD of less than 2 K. 6. Urban and barren surfaces. For UBL, negative biases are obtained for EDLST, positive biases are obtained for the other LST products. All LST products has an RMSD of larger than 2 K. For BRN, except for VLSTO, AGRI-LST, and ABI-L2-LST, the other LST products are overestimated. Only MLST has an RMSD of less than 2 K. 7. All types. EDLST and AGRI-LST are underestimated, especially AGRI-LST; the other LST products are overestimated. The smallest absolute bias is obtained for VLSTO. All LST products have an RMSD of larger than 2 K.
The nighttime comparison demonstrated that the RMSD has a substantial improvement compared with the results during daytime due to the relatively homogeneous distribution in LST during nighttime. The detailed conclusions based on Figure 6b are as follows: 1. Same as the results during daytime, the performance of the MxD11B1 product is highly consistent with that of the MxD11A1 products during nighttime, with an absolute bias of less than 0.5 K and an RMSD of less than 1 K over all surface types. Consequently, we do not further describe the intercomparison results of the MxD11B1 product over each specific surface type. 2. Forests. For ENF, except for AGRI-LST, the other LST products are overestimated. The RMSD of AST_08, AGRI-LST, and ABI-L2-LST is larger than 2 K. For EBF, except for AGRI-LST and ABI-L2-LST, the other LST products has positive biases. The RMSD of AST_08, SL_2_LST, VLSTO, and MLST is less than 2 K; those of the other LST products are larger than 2 K. For DNF, all LST products are overestimated. Only the RMSD of SL_2_LST is larger than 2 K. For DBF, except for AGRI-LST, the LST products have positive biases. Only the RMSD of AGRI-LST and ABI-L2-LST is larger than 2 K. For MXF, EDLST and AGRI-LST have negative biases. Only AGRI-LST and ABI-L2-LST has an RMSD of larger than 2 K. 3. Shrublands and savannas. For CSR, all LST products are overestimated. Except for AST_08, VLSTO, and MLST, the RMSD of the other LST products is larger than 2 K. For OSR, all LST products are overestimated. Only EDLST and SL_2_LST have an RMSD of larger than 2 K. For WDS, except for AGRI-LST, the other LST products are overestimated. The RMSDs of AST_08, AGRI-LST, and ABI-L2-LST are larger than 2 K; those of the other LST products are less than 2 K. For SVN, except for AGRI-LST, the other LST products are overestimated. The RMSDs of EDLST, AGRI-LST, and ABI-L2-LST are larger than 2 K. 4. Grasslands and croplands. For GRS, positive biases are obtained for all LST products. Except for EDLST, MxD21A1, SL_2_LST, and ABI-L2-LST, the RMSD of the other LST products is less than 2 K. For CRP, only AGRI-LST has a negative bias. The RMSD of EDLST, MxD21A1, AGRI-LST, and ABI-L2-LST are larger than 2 K. For CNV, only AGRI-LST is underestimated. Except for AST_08, MxD21A1, AGRI-LST, and MLST, the RMSDs of the other LST products are less than 2 K. 5. Wetlands, snow and ice, and water bodies. For PWL, all LST products are overestimated. The RMSDs of AST_08, EDLST, and MLST are less than 2 K. For PSI, only EDLST has a negative bias. Only the RMSD of AST_08 is less than 2 K. For WTB, only AGRI-LST has a negative bias. The RMSDs of EDLST, SL_2_LST, VLSTO, and ABI-L2-LST are less than 2 K. 6. Urban and barren surfaces. For UBL, all LST products have positive biases. Except for EDLST and MxD21A1, the RMSDs of the other LST products are less than 2 K. For BRN, AGRI-LST and ABI-L2-LST have negative biases. Except for EDLST, SL_2_LST, and AGRI-LST, the other LST products have RMSDs of less than 2 K. 7. All types. Only AGRI-LST has a negative bias. Except for SL_2_LST, AGRI-LST, and ABI-L2-LST, the RMSDs of the other LST products are less than 2 K.
In general, the intercomparison highlighted the key issues related to LST product accuracy assessment: 1. The results for EDLST over most surface types are null because of a large difference of approximately 1 hr between MODIS/Terra (10:30 a.m. and 10:30 p.m. local solar time) and AVHRR/Metop (9:30 a.m. and 9:30 p.m. local solar time) overpass time. 2. Because the ASTER and MODIS sensors are onboarding the same satellite platform (Terra), the LST products of these two sensors could be highly consistent, due to simultaneous observation and nearly the same nadir viewing. In the literature, a consistent performance between ASTER and MODIS LST products has been observed over some homogeneous areas, for example, desert areas (Duan, Li, Cheng, et al., 2017). However, there are still some relatively large disparities between these two products over some land cover types for the intercomparison during daytime as shown in Figure 7. The RMSD over CSR and GRS are higher than 4.6 and 3.2 K, respectively. 3. As shown in Figure 6, compared with MxD11A1, MLST LST is overestimated over all land cover types, especially over ENF, OSR, and WTB. Figure 8 shows the scatterplots of MxD11A1 versus MLST during daytime and nighttime over ENF. The bias and RMSD during daytime are 3.1 and 3.8 K, respectively. Nevertheless, the bias and RMSD during nighttime are significantly smaller than those of the daytime. 4. Many extreme high and low LSTs (e.g., LST >1000 K and LST <70 K) in MxD21A1 are observed over many land cover types. For example, Figure 9 shows the scatterplots of MxD11A1 versus MxD21A1 during daytime and nighttime over BRN. Because no reasonable method is available to eliminate these extreme LSTs, two LST thresholds (i.e., 200 and 350 K) are subjectively selected. Therefore, only the LSTs located between 200 and 350 K in MxD21A1 were used in the intercomparison. The possible reason for these extreme LSTs is the TES algorithm or the inaccurate atmospheric correction. Like MxD21A1, AST_08 and VNP21A1 are also generated by the TES algorithm, but the extreme LSTs are not observed in these two LST products.

Main Problems in Current LST Products
The generation and release of LST products by institutions advance the application of LST products in a variety of studies. To facilitate the use of LST products for broader applications, further research should clarify and solve the problems in LST products. In this section, we identify four main problems in LST products.  1. Spatial discontinuity. As shown in Figure 10, LST values are missing over cloud-covered pixels due to the inability of TIR remote sensing to penetrate clouds, which leads to spatial discontinuity of LST products. 2. Spatiotemporal incomparability. Each pixel in LST products has a different local viewing time and angle. For example, because the MODIS swath width is 2,330 km, the differences in local solar time and VZA for pixels along a given scan line on the same day or for the same pixel on different days in one revisit period could reach 2 hr and 65°, respectively ( Figure 11). Because LST changes with viewing time and angle, LSTs for different pixels on the same day or for the same pixel on different days do not have spatiotemporal comparability. 3. Short time span. As shown in Table 2, except for Landsat and Along Track Scanning Radiometer (ATSR) series LST products with a narrow swath, the single LST products available only date back to the year 2000, which cannot fully fulfill the needs of long time series analysis. However, by combining multiple sources of data to obtain long time series LST product, it is also necessary to consider the significant changes in data observation times due to the orbital drift problem, such as NOAA AVHRR, which make the data incomparable between different times. 4. Instantaneity. Current LST products can only provide instantaneous measurement at the satellite overpass time. However, daily, monthly, and yearly mean LST are required for meteorology-and climatology-related applications. How to estimate daily, monthly, and yearly mean LST from limited instantaneous LSTs is an urgent, meaningful task.
Product users do not concern the advances of the methods for solving these four problems (Sections 5.4.1-5.4.5) are encouraged to skip directly to Section 6.

Filling of LST Over Cloud-Covered Pixels
One of the main deficiencies of TIR remote sensing is its inability to penetrate clouds, which leads to the spatial discontinuity of TIR-based LST products. On average, approximately 68% of the Earth's surface per day is covered by clouds (Stubenrauch et al., 2013), which limits LST applications. To obtain spatially continuous LST products and widen LST applications, accurate filling of LST over cloud-covered pixels is urgently necessary. The various methods that have been proposed to fill LST under cloudy conditions can be grouped into three categories: (a) statistical regression method, (b) spatiotemporal reconstruction method, (c) passive microwave (PMW)-based method, and (d) model-based method.

Statistical Regression Method
The statistical regression method fills cloudy LST in terms of a regression relationship between LST and some auxiliary variables related to LST, for example, the normalized difference vegetation index (NDVI), the digital elevation model (DEM), shortwave radiation, and albedo. The general form of this type of method is written as Various strategies were used to build the regression relationship, for example, multiple linear regression (Y. , artificial neural network (Shwetha & Kumar, 2016), and random forest (Yoo et al., 2020;W. Zhao & Duan, 2020). For determining regression coefficients or the training regression model of Equation 7, two sources of LST were always used, namely, in situ LST and satellite-derived TIR LST. When using in situ LST under cloudy conditions as sample data, the number and spatial representativeness of sample data are crucial for the performance of the statistical regression method. When using satellite-derived TIR LST under clear-sky conditions as sample data, an assumption is that the regression relationship constructed under clear-sky conditions can be extended to cloudy conditions, which could lead to a bias of cloudy LST estimation. To reduce this bias, W. Zhao and Duan (2020) introduced cumulative downward shortwave radiation flux in the regression relationship to characterize the impact of cloud cover on incident solar radiation, which can effectively estimate LST under cloudy conditions. The advantage of the statistical regression method is simple, and its disadvantage is area-dependent.

Spatiotemporal Reconstruction Method
The main idea of the spatiotemporal reconstruction method is that cloudy LST is considered the sum of theoretical clear-sky LST and a correction term of LST. This method is expressed as where LST cld is the filled LST under cloudy conditions, LST clr is the theoretical clear-sky LST, and ΔLST is the LST correction term.
The spatiotemporal reconstruction method comprises two steps. The first step is the reconstruction of clear-sky LST. In the literature, many methods have been developed to reconstruct clear-sky LST in terms of information from adjacent pixels, including temporal information (Y. Xu & Shen, 2013), spatial information (Neteler, 2010;, and spatiotemporal information (Hengl et al., 2012;Metz et al., 2014;L. Sun et al., 2017;Yao, Wang, Huang, Sun, et al., 2021). These methods can effectively reconstruct clear-sky LST. However, the reconstructed LST is theoretical clear-sky LST, which differs from actual cloudy LST.
The second step of the spatiotemporal reconstruction method is the estimation of the LST correction term. According to different parameterization schemes, two methods were developed to estimate the LST correction term, namely, the net shortwave radiation difference method (Jin, 2000;Jin & Dickinson, 2000) and the air temperature difference method (Z. Zou et al., 2018). For the net shortwave radiation difference method, the correction term of LST was parameterized as a function of the difference in net shortwave radiation between a cloudy pixel and its adjacent clear-sky pixels (Jin, 2000). For the air temperature difference method, the correction term of LST was parameterized as a function of the difference between actual air temperature and theoretical air temperature modeled by an annual temperature cycle model (Zou et al., 2018). These two methods still require other auxiliary variables, for example, wind speed, SM, and relative humidity, from weather station observations or reanalysis data as input data, which probably limit their applications over extended areas.

PMW-Based Method
The main idea of the PMW-based method is merging clear-sky TIR LST and cloudy LST from PMW data for the generation of all-weather LST products. TIR measurements are widely used to retrieve LST with high accuracy and high spatial resolution but are limited to clear-sky conditions due to their inability to penetrate clouds. By contrast, PMW measurements can penetrate clouds and provide data regardless of the cloud conditions. However, PMW measurements have their limitations, such as low spatial resolution and low temperature-retrieval accuracy . Furthermore, temperature retrieved from PMW measurements is the subsurface temperature with a penetration depth approximately ranging from 1 to 10 mm, depending on the wavelength used to measure, which differs from the skin temperature retrieved from TIR measurements with a penetration depth approximately ranging from 1 to 100 μm (Hecht, 2002).
Because different characteristics and limitations make TIR-and PMW-derived LST complementary, the idea of combining these two types of LST products to generate an all-weather LST product has been proposed (Aires & Prigent, 2006;Holmes et al., 2015; Z.-L. Li et al., 2013;Parinussa et al., 2008). For example, Kou et al. (2016) attempted to merge TIR-and PMW-derived LST by using the Bayesian maximum entropy method. A framework for the generation of an all-weather LST product by merging TIR-and PMW-derived LST was developed by . Since then, studies have attempted to merge TIR-and PMW-derived LST using various fusion strategies, for example, temporal component decomposition (X. , cumulative distribution function matching, and multiresolution Kalman filtering (S. Xu & Cheng, 2021). Figure 12 shows a framework for the generation of an all-weather LST product in combination with TIR and PMW data. In this process, two key steps are involved: (a) conversion between skin temperature and subsurface temperature, and (b) spatial downscaling of cloudy PMW-derived LST. The two steps have influences on the accuracy of the all-weather LST product. Although some studies have attempted to convert between skin temperature and subsurface temperature F.-C. Zhou et al., 2016;J. Zhou et al., 2017), this problem has not been well solved, because no physical-based model effectively describes variations in the soil temperature profile. To merge TIR-and PMW-derived LST, it is necessary to spatially downscale cloudy PMW-derived LST for matching the spatial resolution of clear-sky TIR LST. Discovering how to effectively characterize spatial variations in cloudy LST within a PMW pixel by using auxiliary parameters at high spatial resolution is the key to solving this problem. In addition to the two aforementioned steps, the retrieval accuracy of PMW-derived LST has a significant impact on the accuracy of the all-weather LST product. Duan et al. (2020) conducted an exhaustive review on the methods of LST retrieval from PMW measurements. They demonstrated that the retrieval accuracy of PMW-derived LST is limited by the estimation accuracy of PMW surface emissivity because PMW emissivity varies strongly with surface properties, such as SM, vegetation cover, and surface roughness.

Model-Based Method
With the development of land surface modeling and data assimilation technique, several land surface models (LSMs) in well-established land data assimilation system such as the Global Land Data Assimilation System (GLDAS) and North American Land Data Assimilation System (NLDAS) can output spatially complete and temporally continuous LST. In addition, reanalysis data set (e.g., ERA5-Land) from numerical weather prediction (NWP) models also provide spatially complete and temporally continuous LST. Similar as the PMW-based method, the model-based method also merges TIR In general, the aforementioned methods are promising for reconstruction of cloudy LST. Because of the complementary characteristics of TIR-and PMW-derived LST, combining the advantages of these two LSTs is attractive for generating an all-weather LST product at TIR spatial resolution, especially over large cloud-covered areas.
Although several methods were proposed in this aspect, many substantive problems have not been solved; thus, improvements could be made. First, the retrieval accuracy of PMW LST should be improved. This endeavor requires higher calibration accuracy of PMW sensors and the development of new physically based retrieval algorithms of PMW-derived LST, for example, the SW-like PMW LST retrieval algorithm or the simultaneous retrieval algorithm of PMW-derived LST and SM. Second, spatial downscaling of PMW-derived LST under cloudy conditions should be refined. Unlike some traditional methods, machine learning-based methods (e.g., random forest and lightGBM) could show a promising perspective in the spatial downscaling of cloudy PMW-derived LST. Last, developing a physical-based method is necessary for the conversion between subsurface temperature and skin temperature in terms of a soil heat conduction equation.

Temporal Normalization of LST
As shown in Figure 11, there are obvious differences in local solar time for pixels in the east-west direction along the scanning line because of the large swath width in many satellite observations. Taking the MODIS sensor as an example, because of its viewing swath width of 2,330 km, there are more than 1 hr difference in local solar time at the west-east direction. This scanning feature results in unavoided temporal inconsistency in current LST estimates from wide-field-of-view scanners. To increase the spatiotemporal comparability of LSTs for different pixels on the same day or for the same pixel on different days, many efforts have been made to correct for the temperature difference induced by the temporal difference through normalizing the LST estimates to the same local solar time. According to the temporal distribution of the clouds within a day, the normalization methods can be divided into two typical categories: the diurnal temperature cycle (DTC)-based normalization method for a fully cloud-free condition and the empirical normalization method for a partial cloud-free condition.

DTC-Based Normalization Method
The diurnal LST temporal evolution on a cloud-free day can be described by a DTC model. Based on a different theoretical basis, two typical DTC models have been developed: the thermal inertia (TI)-based DTC model and the semi-empirical (SE)-based DTC model.
• TI-based DTC model LST dynamics are coupled into surface energy balance (SEB) process through energy exchanging terms such as ground flux, sensible heat, and latent heat. According to near-surface conductive heat transfer, the temperature variation obeys heat conduction equation. To solve this equation and depict LST evolution within a day, many studies have used specified boundary conditions, namely, surface energy balance and constant soil temperature at a certain depth. Through solving the equation with the Fourier series method and using different simplification schemes for SEB, the general TI-based DTC model was mathematically expressed as follows (Price, 1977;Sobrino & El Kharraz, 1999;Watson, 1975;Xue & Cracknell, 1995): where T s (t) is the LST at time t, and h 0 , h 1 , and thermal inertia (P) are the unknown (free) parameters. Other parameters (A n , B n , ω, and δ n ) can be assessed with solar angle, geolocation, and date information. To determine the three unknowns, there should be no less than three LST measurements. Because the first term (n = 1) dominates the temperature series (Price, 1977), this TI-based DTC model is usually simplified as the first-order expression (Xue & Cracknell, 1995). However, the experiment results indicated that it prefers using at least a two-order (n ≥ 2) expression during daytime and a three-order expression (n ≥ 3) at night to obtain a high-accuracy description of the diurnal LST changes.
In theory, the TI-based DTC model can describe the diurnal temperature evolution by using only three LST measurements within a day. The limitations of this model include: (a) the high complexity of the model induces its numeric solution unstable with limited LST measurements; and (b) the initial values of the free parameters h 0 and h 1 are difficult to assign because of their unclear physical meanings.
To alleviate the limitations of the aforementioned TI-based DTC model, many SE-based DTC models have been proposed. On the basis of the thermal diffusion equation (heat conduction equation), the common structure of SE-based models is a harmonic function to predict the evolution of the daytime LST (T day ) in conjunction with an exponential, hyperbolic, polynomial, or harmonic function to describe the decay of the LST at night (T night ). Additionally, the daytime and nighttime LST functions should be continuous and derivable at the transition time of the two functions. A typical SE-based DTC model was expressed as follows (Göttsche & Olesen, 2001): The meanings of all the parameters in the model are illustrated in Figure 13, where T 0 is the residual temperature around sunrise; T a is the temperature amplitude; t m is the time of the maximum temperature; t s is the start time of free attenuation; δT is the temperature difference between T 0 and T (t→∞); k is the attenuation constant of LST at night; and ω is the width over the half-period of the cosine term, which can be determined by the solar geometry or the simplification of the heat conduction equation at time t sr (Duan et al., 2013). Therefore, to derive the five unknown (free) parameters (T 0 , T a , t m , t s , and δT), there should be no less than five LST observations available within a day.
Subsequently, many other SE-based DTC models have been proposed (Göttsche & Olesen, 2009;Inamdar et al., 2008 Zhan et al., 2012) that focus on fully filling the demand for accurate fitting of the LST diurnal cycle. These SE-based DTC models usually have five or six unknown parameters. Nevertheless, for commonly used polar-orbit satellite observations, such as the MODIS LST product, there are usually no more than four observations per day at a given location (except for the polar regions), and the free parameters in the SE-based DTC model cannot be well solved. To limit the degrees of freedom and increase conversion to a solution, many different parameter-reduction approaches have been proposed to obtain a priori knowledge of at least one or two parameters. For example, t s was approximated as a fixed value of 1 hr before the time of sunset (Duan et al., 2014b), Holmes et al. (2013) expressed t s as a function of t m , ω, T a , and δT,  set the optimal value of the total atmospheric optical thickness in the SE-based DTC model developed by Göttsche and Olesen (2009) to 0.01. In addition to reducing the number of unknown parameters, the incorporation of multi-source remote sensing data is used to assist the DTC fitting by increasing the number of observations or by providing a priori knowledge of some model parameters. For instance, the geostationary observations could provide the information of t m (D. Sun & Pinker, 2005;J. Zhou et al., 2011).
In general, the aforementioned two types of DTC models provide different options for depicting diurnal LST temporal evolution in a day under a cloud-free condition. By selecting the right model and effectively fitting the model with limited LST observations in 1 day, the model's free parameters can be determined in a fully cloud-free condition. Once the DTC model determined, it can be used to describe the diurnal LST temporal evolution, obtain LST at any given local solar time, and realize the LST temporal normalization.

Empirical Normalization Method
Different from the DTC-based normalization method with its requirement of an accurate description of temperature variation at a diurnal scale, the empirical regression method mainly uses the temperature variation pattern at a certain period in a day to perform temporal normalization. As shown in Figure 13, the LST changes in the mid-morning (10:00-12:00) approach a linear increasing process in a cloud-free condition. According to this Figure 13. Illustration of the SE-based DTC model parameters. T day (t) and T night (t) are LSTs of the daytime and nighttime parts, respectively. t is the time, T 0 is the residual temperature around sunrise, T a is the temperature amplitude, ω is the width over the half-period of the cosine term, t m is the time at which the temperature reaches its maximum, t sr is the time of sunrise, t s is the start time of free attenuation, δT is the temperature difference between T 0 and T (t→∞), and k is the attenuation constant (Adapted from Duan et al. [2012]). feature, Duan et al. (2014a) proposed an empirical regression method to correct the temperature differences induced by the time difference based on the LST changing slope (T slp ) within this period: where T s (t) is the satellite-derived LST at time t, t N is the normalized time, and T s (t N ) is the temporally normalized LST. Duan et al. (2014a) expressed T slp as a function of NDVI, the cosine of solar zenith angle (SZA), and surface elevation: where d 0 -d 3 are the coefficients to be determined in advance, θ s is SZA, and Elv is the surface elevation. Figure 14 shows an example of the LST temporal normalization with this method. Obviously, the original LST product has a clear temperature difference (Figure 14a) induced by the difference in local solar time (Figure 14d). After normalization, the spatial variations in the LST (Figure 14b) appear to be more uniform than those in the original product. As shown in Figure 14c, the maximum corrected differences among different pixels can reach to 6 K. In addition to the aforementioned derivation of the changing slope, there would be other options such as the estimation with the use of the slope estimated from the BT of geostationary satellite observations (Duan et al., 2013).
It should be noted that the method proposed by Duan et al. (2014a) is only appropriate for conducting LST temporal normalization in the mid-morning period. For LST observations at approximately noon (i.e., Aqua-MODIS), the linear approximation is not suitable. However, the quadratic polynomial would be a satisfactory choice to represent the temperature variation during this period. Through constructing a correction coefficient via the quadratic polynomial function, the original LST observations from different local solar times can be normalized to the values at the reference time.
In addition to numerically fitting the LST variations at a certain period (middle-morning or noon) with the assistance of other satellite observations, the introduction of the LST simulations from land surface model (LSM) is also an alternative. Because they are driven by local atmospheric forcings, LSM could provide LST diurnal pattern at different land cover conditions and the simulations should have similar LST temporal changing features as those of the pixels in real circumstances. By fitting the changes in the LST simulations for each pixel (Jin & Dickinson, 1999), the temporal pattern derived from the simulated LSTs can be transformed to derive the actual LST changing feature of each pixel. Then, this temporal information of LST changes can be applied to normalize the LST observations from satellites to obtain the LSTs at a consistent time. In general, with the aid of LSM simulation, these short-term variation patterns in the LST diurnal cycle provide high flexibility to temporal normalization.
Although the aforementioned empirical normalization approaches have used numerical expression or the LSM model to express the relationship between LST changes and surface factors (vegetation, topography, and solar illumination), the strong coupling effect of LST with other surface factors makes the relationship expression complicated, even for the LST simulation. Comparatively, the machine learning method has the advantage of modeling the complex non-linear relationship. On this basis, W.  proposed a machine learning method to build a linking model to express the complex relation between LST and the solar radiation factor and other surface factors such as the vegetation index, topographic factors, and the wetness index. Under the prerequisite that the temporal inconsistency in LST observations mainly resulted from the differences in solar radiation induced by the time difference of satellite observation, the temporally normalized LST was then obtained by driving the linking model with a temporally consistent solar radiation factor. This method is applicable to any satellite observation during the daytime but requires knowledge of the solar radiation factor in advance.
In conclusion, many efforts have made positive effects on deriving temporal consistent LSTs. The incorporation of multi-temporal observations in a day from geostationary satellite or constellation offers important information to describe LST diurnal variation and normalize the temporal differences. However, the DTC-based normalization method is mainly conducted under the assumption of clear-sky conditions, which is not satisfied in most circumstances. The appearance of cloud cover during the daytime strongly distorts the LST temporal evolution and makes the predefined DTC models useless. Although the empirical normalization method at the middle-morning period slightly lightened this impact by using the clear-sky assumption within a short period, the accuracy of the coefficients in Equation 13 directly influences the normalization accuracy. Comparatively, the physical modeling of LST evolution at a certain period with LSM provides a satisfactory option for temporal normalization. However, the normalization is not sufficiently accurate and is highly dependent on the accuracy of the model simulation. The machine learning-based empirical normalization does not have the clear-sky assumption and can partly resolve the complex interaction between LST and other surface parameters. However, the correction suffers from the accuracy of solar radiation and shows high dependency on the training data according to the characteristics of the machine learning method. In summary, although current studies have made attempts to address this issue, there are still plenty of spaces in the fields of LST temporal normalization. How to develop a more physical model to depict LST diurnal evolution in all-weather conditions but with limited degrees of freedom is a substantial challenge for further research in this field.

Angular Normalization of LST
In natural conditions, the land surface is always heterogeneous at a satellite pixel, which consists of several homogeneous components with different temperatures (e.g., sunlit and shaded soil/vegetation). As shown in Figure 15a, these components have different area weights under different solar and observation geometries, leading to different observed directional LSTs. Various experiments, including those at ground (Z.-L. Li et al., 2004), airborne (Lagouarde et al., 2010), and satellite scales , reported this angular anisotropy and revealed that the LST variation can reach up to 15 K ( Figure 15b). Evidently, the angular anisotropy of LST would degrade the comparability of LST and ultimately restrict wide applications of LST (Rasmussen et al., 2011). Therefore, it is strongly necessary to conduct angular effect correction on the remotely sensed LST, namely, normalizing LSTs observed in different angles to the value in a reference angle, for example, the nadir direction. Some explorations have been conducted on the angular normalization of LST, and two types of methods have been proposed: the kernel-driven model-based method and the component temperature aggregation method.

Kernel-Driven Model-Based Method
The kernel-driven model assumes that the LST angular effect is a linear combination of different kernels, which can be formulated as follows: where AE LST is the angular effect of LST, which can be expressed by LST (X. Liu et al., 2021), and the LST difference (Duffour et al., 2016) or LST ratio (Vinnikov et al., 2012) between the observed direction and reference direction; 1, K geo , and K vol are the isotropic kernel, geometric-optical (GO) kernel, and volume scattering kernel, respectively; these kernels are mathematical functions of observation angles, solar angles, and canopy structural parameters, representing component fractions or fraction differences (Wanner et al., 1995); f iso , f geo , and f vol are kernel coefficients, corresponding to different component temperature or temperature differences (X. Liu, Tang, Li, Zhang, & Shang, 2020); θ v , θ s , and φ are VZA, SZA, and relative azimuth angle, respectively.
The current kernel-driven models in the TIR domain are summarized in Table 3. Notably, the kernel-driven model was originally designed for simulating the bidirectional reflectance distribution function (BRDF) in the visible and near-infrared domain. Because reflectance and thermal radiation share the same geometric structure, Lagouarde & Irvine (2008) adapted an optical hotspot model (Roujean, 2000) to the TIR domain by replacing reflectance with temperature, and developed RL model. Similarly, Ren et al. (2014Ren et al. ( , 2015 extended a classic BRDF model, that is, Ross-Li model (Wanner et al., 1995), to the TIR domain by replacing reflectance with thermal radiance. The main limitation of these two BRDF-developed models is that the RL model does not have K vol , and the K vol of Ross-Li model is not suitable for TIR radiation . By contrast, L. Su et al. (2002) and Vinnikov et al. (2012) have respectively proposed a kernel-driven model that describes the influences of emissivity directionality. The former proposed the LSF-Li model derived from a conceptual model for effective directional emissivity, and the latter proposed the Vinnikov model, which is based on a comprehensively statistical analysis of synchronous LST products from two geostationary satellites. In addition, unlike BRDF modeling, based on temperature differences rather than illumination differences (Wanner et al., 1995), X. Liu, Tang, Li, Zhang, and Shang (2020) established a novel K geo that assumes that the sunlit soil component has the highest temperature and that other component temperatures are the same and then proposed a GOg-Emi model and a GOg-LSF model by considering the directional emissivity effect (X. Liu et al., 2021). The aforementioned three-parameter models either have poor overall performance (RL model) (Cao et al., 2019;X. Liu et al., 2018) or would underestimate the hotspot effect (Ross-Li, LSF-Li, Vinnikov, GOg-Emi, and GOg-LSF model) (Cao et al., 2019;. Correspondingly, some four-parameter models have been developed by combing the advantages of different kernel functions Ermida, Trigo, DaCamara, & Roujean, 2018). Their main characteristic is using an extra unknown variable (hotspot width) to improve the fitting capability in the hotspot region. Figure 16 illustrates the flowchart of the angular normalization of LST based on kernel-driven models. First, developing an LST multi-angular data set with adequate angle samples to solve kernel coefficients is necessary. Once the model is calibrated, namely, the free parameters are determined, LST can be obtained at reference direction with the calibrated model, thus completing the angular normalization of LST. However, current kernel-driven models are three parameters or four parameters, whereas existing products (  (2018) used one polar-orbit satellite and one geostationary satellite, and X. Liu et al. (2018) combined two polar-orbit satellites. However, even with these efforts, only the dual-angle LSTs of the same pixel can be obtained, and the solution of the kernel-driven model remains not obtained. As a result, some hypotheses that weaken spatial and temporal heterogeneity have been proposed to enrich angle information, for example, that LSTs of the same pixel have the same angle effect in the same year (Ermida et al., 2017) and that LSTs in the same land cover type (X. Liu et al., 2018)   have the same angular effect. Figure 17 is an example of the angular normalization of LSTs and shows that the spatial variations in the LST after angular correction (Figure 17b) are more uniform than those in the LST before angular correction (Figure 17a). In addition, LST differences ( Figure 17c) are highly related to VZA differences ( Figure 17d). However, the validity of the aforementioned assumptions on the consistency of LST angular effects can exert an important influence on the eventual normalization performance. In addition, LSTs from different satellites have a systematic error, and this error may be larger than the LST angular effect.
In summary, the problems with the kernel-driven-model-based method are the absence of an operational but high-accuracy model and the limited multi-angle information used for model calibration. In the case of kernel-driven modeling, more TIR characteristics that are distinct from the visible and near-infrared, for example, temperature inertia effect (X. , should be considered. In addition, improving the model performance (especially in hotspot regions) without adding extra parameters remains a problem to be solved. Regarding model calibration, the use of more (i.e., ≥3) polar-orbit satellites with a large scan swath may provide a more complete angle coverage. Moreover, manufacturing a new TIR sensor with the multi-angular configuration designed by the theoretical analysis of kernel-driven models (Ran et al., 2021;Ren et al., 2015) is a promising method.

Component Temperature Aggregation Method
Temperature contrast among components is the main driver of the angular effect of LST. Therefore, if the temperatures of different components in the mixed pixel can be separated from the retrieved LST, the LST at the reference direction can be estimated by aggregating component thermal radiances to pixel thermal radiance according to a proper thermal radiation directionality (TRD) model. The flow chart of the component temperature aggregation method is shown in Figure 18. First, different available information is used to separate component temperatures with the aid of a TRD model. This TRD model should be a retrievable, accurate analytical model. Next, component temperatures are aggregated to pixel LST at the reference direction based on a TRD model. This model can be the aforementioned TRD model for separation, or a more complex (i.e., no explicit expressions) but comprehensive physical model. Therefore, auxiliary data such as canopy structure parameters may be required.
Relevant studies have mainly explored methods for separating component temperatures, especially for two components (i.e., soil and vegetation). According to available data sources, current efforts can be classified into the multi-angle method (Jia et al., 2003;Z.-L. Li et al., 2001), multi-channel method (Sobrino & Caselles, 1991;Xie et al., 2016), and multi-pixel method (Song et al., 2015;W. Zhan et al., 2013;R. Zhang et al., 2005), corresponding to the use of angular, spectral, and spatial information, respectively. Although the multi-angular method is an effective method for the separation of component temperatures, the absence of multi-angular TIR images limits its wide applications. Until now, only the ATSR series and Sentinel-3/Sea and Land Surface Temperature Radiometer (SLSTR) provided two angular TIR observations at nadir and forward 55°. The multi-channel method can overcome the problem of data deficiency, but the high correlation between different channels extremely restricts the utility of this method. By contrast, the multi-pixel method has satisfactory feasibility. However, it assumes that component temperatures within a neighborhood (e.g., 5 × 5 pixels) are consistent, which may increase the difficulty to obtain a computationally stable solution when correlations between pixels are high (X. Liu, Tang, Li, Zhou, et al., 2020;W. Zhan et al., 2013). Because of the limitation of using a single information source, some improved methods that combine multiple available data sources have been proposed. Although the availability of separation methods continues to be improved, the estimation error was still remarkable, larger than 2 K (X. Liu, Tang, Li, Zhou, et al., 2020), which cannot support the accurate angular correction of LST.
Regarding angular normalization, there are a few exploratory studies. For example, X. Liu, Tang, Wu, et al. (2019) retrieved soil and vegetation component temperatures using the multi-pixel method calculated component proportions with the help of the BRDF model and finally constituted angular normalization of LST according to a GO model. Analyses based on a simulated data set revealed that angular normalization can improve the LST retrieval accuracy caused by angular effect, from 1.2 to 0.8 K. The main limitation of this method is that only two components are considered. Other temperature differences, especially the temperature contrast between sunlit and shaded components, are ignored, significantly degrading the correction ability.
The aforementioned discussion indicates that the problem with the component temperature aggregation method is the absence of a feasible and accurate TRD model. Regarding component temperature separation, the TRD model should be sufficiently retrievable and accurate. In addition, it could be used to estimate temperatures of three (i.e., vegetation and sunlit/shaded soil) and even four components (i.e., sunlit/shaded soil and sunlit/ shaded vegetation). At this point, the GO model (Pinheiro, Privette, & Guillevic, 2006) with some reasonable parameter-reduction scheme (X. Liu, Tang, Li, Zhang, & Shang, 2020;Wanner et al., 1995) is a possible method. Regarding component temperature aggregation, the desired TRD model should be sufficiently comprehensive to describe the relationship between the component temperature and pixel LST accurately. Moreover, other variables can be obtained in practice to ensure the applicability of this model.

Long-Term Time Series of the LST Product
Many application fields, such as climate change, urgently need over 20-30 years of globally covered and high spatiotemporal resolution (e.g., 1 km daily) LST products (Ghent, 2014). Although the Landsat and GOES series of satellites can provide long-term historical LST image, the temporal resolution of Landsat data is only 16 days, and GOES, as a geostationary meteorological satellite, cannot achieve global coverage, and early data is difficult to obtain, and neither of them can fully fulfill the demand. However, current LST remote sensing products with global coverage and daily observations only date back to the year 2000.
Fortunately, NOAA AVHRR data, which have recorded easily accessible historical global top-of-atmosphere BT twice daily in two TIR channels with 1.1 km spatial resolution from 1981 starting with the NOAA-7 satellite (Pinheiro, Mahoney, et al., 2006;Z. Qin, Dall'Olmo, et al., 2001), are currently treated as the most frequently used data source to generate long-term LST datasets ( Nonetheless, owing to the lack of active controllers that remain on their nominal sun-synchronous orbits, the crossing time of NOAA satellites shifts backward gradually after their positioning, which is on average approximately 30 min per year (Price, 1991). The pattern of Equatorial Crossing Time of each NOAA satellite is shown by the solid lines in Figure 20, and the longest change is more than 3.5 hr. Orbital drift leads to the variability for multi-temporal TIR remote sensing images in the local solar time and viewing angles (Pinheiro et al., 2004), making the LST correction induced by orbital drift a necessary pre-processing step for long-term LST analysis when using AVHRR data. Two categories of correction methods, namely, the statistical regression method and the LST daily cycle reconstruction method, have been proposed to correct for the orbital drift effect.
The statistical regression method selects a proxy signal, usually, the SZA, to represent the image acquisition time (Julien and Sobrino, 2021) and then fits an empirical function between the proxy signal and the correction values of the LST from the time series data (Gleason et al., 2002;Gutman, 1999;Julien & Sobrino, 2012;Sobrino et al., 2008). An example of the fitting function is shown in Equation 15: where Δ LST (t) is the LST correction value, Δ SZA (t) is the SZA difference before and after correction, t is the observation time, and α i (i = 0, 1, 2) is the fitting coefficients. The statistical method is computationally simple, has few input parameters, and can be directly used for large-scale images. However, the correction results are empirical estimates without physical meaning.
By contrast, LST daily cycle reconstruction corrects the LST of each pixel for different observations to the same reference time by building the DTC model described in Section 5.4.3.1, which is a physical method suitable for where Δ LST (t) is the LST correction value, T DTC is the LST described using the DTC model, t is the original observation time, and t u is the uniform observation time after correction. Thus, the correction results are physically meaningful, but their performance is highly related to the accuracy of the DTC model, which requires multiple observations.
After orbital drift correction, the AVHRR LST product can be combined with the MODIS LST product to generate a long time series (more than 40 years) LST product. However, as an alternative to compensate for the lack of data, the root mean square error (RMSE) of the AVHRR LST after orbital drift correction is approximately 3 K according to the validation results (Julien & Sobrino, 2021;X. Liu, Tang, Yan, et al., 2019), which is less accurate than the LST product generated by the original TIR remote sensing data. In addition, for the long-term time-series of LST products constructed with multiple data sources, there are high requirements for the calibration accuracy of individual sensors and consistency between different data sources to ensure the data quality and application reliability (Xiong et al., 2020). For sensors that have been in operation for a long time, such as MODIS, corrections are needed to eliminate the effects of instrument degradation (Lyapustin et al., 2014), while for combinations of different sources, the use of inter-calibration is considered to reduce systematic errors caused by the diversity between sensors (Duan, Li, Cheng, et al., 2017;S. Li & Jiang., 2018;Ye et al., 2021).

Temporal Aggregation
Existing LST products typically only provide clear-sky temperature at satellite overpass time. Some institutes further average these instantaneous temperatures to a monthly scale for every nominal overpass time (e.g., MxD11C3 of the National Aeronautics and Space Administration). However, a diversity of studies, such as those on climate change and meteorology, actually need true mean LST at a daily/monthly/annual scale (i.e., an average of all instantaneous LSTs regardless of whether they are cloudless). Therefore, temporal upscaling of instantaneous LST to mean LSTs for a specific period can effectively bridge the gap between existing products and actual requirements.

Daily Mean LST
Nowadays, algorithms for deriving daily mean LST (DMLST) can be divided into three major categories ( Table 4).
The first method is based on the DTC models. As aforementioned, DTC models can simulate instantaneous LST for any moment; thus, obtaining DMLST is easy as long as the free parameters in DTC models have been determined. However, the problem is how to select suitable models and solve frameworks for polar satellites due to their sparse temporal samples. If four measurements are available, several four-parameter DTC models (Duan et al., 2014b;Hong et al., 2018) can be used directly to obtain DMLSTs. However, if this ideal situation cannot be fulfilled or a DTC model with more than four parameters is used, the incorporation of other information, for example, neighborhood observations, component temperatures (X. Liu, Tang, Wu, et al., 2019;Quan, Chen, Zhan, Wang, Voogt, Li, et al., 2014), and intra-annual temperature variation (Hong et al., 2021), can make ill-posed models solvable. Notably, DTC models are derived from all clear-sky conditions, which leads to two shortcomings of DTC-based DMLST: the total number of DMLSTs results is relatively small because all overpass moments being cloudless is rare in reality, and estimated DMLST has a positive bias relative to the true DMLST, which is approximately 1 K according to the findings of Hong et al. (2021).
The second and third methods are based on statistical regression models and directly use true DMLST as the target variable. These two methods learn from the calculation of the mean air temperature. When temporal-continuous measurements are not available, which is similar to the sparse observations of polar satellites, a common method is estimating daily mean air temperature by using a linear combination of data collected at several special moments (Nordli et al., 1997). In the context of DMLSTs, a widely used model is the simple average of two instantaneous LSTs (denoted as TSA) of afternoon satellites, for example, MODIS (Alkama & Cescatti, 2016;Compo et al., 2013;Duveiller et al., 2018;Y. Li et al., 2016; and Atmospheric Infrared Sounder (Ding et al., 2020;Susskind et al., 2019). The basis is that the overpass times of afternoon satellites (i.e., 13:30 + 01:30) can be considered the times of maximum and minimum LST (Beale et al., 2020). However, there are differences between actual overpass time and nominal overpass time. The times of maximum and minimum LSTs also vary due to many factors, such as cloud interference and SM changes. In addition, the simple average of two clear-sky LSTs may generate a more positive bias than true DMLST because the clear-sky nighttime LST is probably not the minimum LST. These limitations may reduce the validity of the TSA model. Correspondingly, Xing et al. (2021) proposed to use the weighted average of multiple instantaneous LSTs (denoted as MWA) to estimate DMLST. This model includes nine combinations in which there is at least one daytime moment and one nighttime moment, and their time difference is approximately 12 hr. To calibrate this practical model, Xing et al. (2021) collected numerous field measurements from 235 stations worldwide and then developed the linear relationship between true DMLSTs including cloud and cloudless conditions and clear-sky instantaneous LSTs of satellite overpass time. Validation results reveal a higher accuracy of the MWA model (RMSE less than 1.6 K) than the TSA model (RMSE ∼ 2.6 K). Figure 21 shows an example of global DMLST distribution estimated by the TSA model and the MWA model, and their differences. TSA-based DMLST is usually higher than MWA-based DMLST (Figure 21c). In addition, comparing the subplots in Figures 21a and 21b shows that a larger number of moment combinations for the MWA model also produced a larger number of valid DMLSTs than the TSA model.  (2016) MWA model-based The result is true DMLST, with high accuracy and complete spatial coverage The universality of the model is significantly affected by the richness of samples Xing et al. (2021) Note. TSA is a simple averaging of two observations, and MWA is a weighted averaging of multiple observations.  In summary, the MWA model has satisfactory operability, high accuracy, and complete spatial coverage and therefore is a promising method for estimating DMLST. However, its performance largely depends on the validity of weighted coefficients of different times. The uneven distribution of ground stations used to develop the model, especially the data deficient in Africa and South America, may degrade the robustness of the MWA model. Therefore, collecting comprehensive in situ measurements or combining the continuous observations of geostationary satellites to obtain reliable model coefficients is necessary.

Monthly Mean LST
Similar to DMLST, there are three models for conversing instantaneous LST to monthly mean LST (MMLST), namely, the DTC model, TSA model, and MWA model. To calculate MMLST, the DTC model is applied to monthly data. Correspondingly, this extended version of the DTC model is called the monthly mean diurnal cycles (MDC) model (Ignatov & Gutman, 1999). Specifically, observations on different days within the same month are first averaged according to overpass time and then used to develop one "DTC" model for this month. For this reason, MDC model requires that the number of average LSTs in different moments is no less than four, rather than every instantaneous LST of every day being under the clear sky. Similar to DTC models, there are two situations for MDC models: use a four-parameter model or appropriate interpolation method (Aires et al., 2004;Sharifnezhadazizi et al., 2019) and synthesize more useful information to compensate for the mismatch between the number of observation averages and the number of model coefficients, such as fusion observations from multiple instruments (Holmes et al., 2015;Norouzi et al., 2015) and a priori knowledge retrieved from geostationary satellites . Regarding the TSA model, two instantaneous LSTs from morning satellites, for example, the ATSR series (Good et al., 2017) and MODIS aboard on Terra (X. L. , have also been applied to retrieve MMLST, when afternoon observations are unavailable. For the MWA model, the situation is the same as that of DMLST estimation. Different from DMLST, for MMLST estimation, further selecting a temporal aggregation scheme is necessary. The current methods, the DTC model-based (i.e., MDC) and the statistical model-based (i.e., TSA and MWA), first obtain the LST averages at each observation time and then develop the related models using these averages. This scheme is called average by observations (ABO). Another approach is called average by DMLSTs (ABD). Specifically, estimate DMLSTs and then average all available DMLSTs to calculate MMLST. The schematic for these two schemes is illustrated in Figure 22. Obviously, the main difference between these two schemes is that the ABO approach assigns equal importance to all participating instantaneous observations within a month, whereas ABD treats all DMLSTs equally (Ding et al., 2020). Figure 23a is the comparative result of six methods for estimating MMLST, which are the cross-combinations of the aforementioned two temporal aggregation schemes and three conversion models by using the same data as Xing et al. (2021). This result reveals that the influence of the conversion model is predominant (RMSE differences larger than 0.7 K), whereas the effect of the temporal aggregation scheme is slight (RMSE differences less than 0.2 K). Furthermore, for the MWA model, although the overall accuracies are the same for using the ABO and ABD scheme, the combination with ABD can result in a better performance than that with ABO when the number of DMLST (NOD) is ≥20 (see subplot b). Therefore, the best strategy for estimating MMLST can be determined using the MWA model and ABO when the NOD is <20, whereas using MWA model and ABD when NOD is ≥20, and the RMSE is approximately 1.5 K. Figure 24 shows the distribution of global MMLSTs in 2018-07 using this combination strategy and MOD11A1/ MYD11A1 products. When comparing Figures 21b and 24, it can be observed that MMLSTs have more complete spatial coverage than DMLSTs. Specifically, there are little data missing in the rainforest area. This result indicates that the temporal upscaling can fill in the spatial gap effectively. However, similar to the performance of DMLST, the performance of the best strategy for MMLST is also significantly influenced by the robustness of the MWA model.

Annual Mean LST
Regarding the annual mean LST (AMLST), a prevalent method is using annual temperature cycle (ATC) models because AMLST is a main parameter within the ATC models. Different from DTC and MDC models, ATC models only use sine or cosine functions to fit the LST annual variation. Two representative models are the ATC model with three free parameters (ACP3), namely, Equation 17 (Bechtel, 2011), and the ATC model with five free parameters (ACP5), namely, Equation 18 (Bechtel & Sismanidis, 2018), such that where doy is the day of the year; amp 1 is the annual amplitude of DMLST; ph 1 is the phase, namely, the day when DMLST reaches its yearly maximum; ω is the annual period, which is 365 (366) for a common year (leap year); amp 2 and ph 2 are the amplitude and phase of biannual variations for DMLSTs, respectively. The inputs of ATC model are DMLSTs on different dates within the same year, and the model free parameters are usually three (ACP3 model) or five (ACP5 model). Therefore, the number of observations is no longer an obstacle for the ATC model. By contrast, selecting an appropriate model is the most important. For instance, based on rich sample points selected every 1° latitude from the 27°E line, Xing et al. (2020) compared the abilities of the ACP3 model and the ACP5 model (Bechtel & Sismanidis, 2018) to describe annual LSTs. The result is shown in Figure 25, indicating that the ACP5 model has an obviously better performance than the ACP3 model only in tropical and Antarctic, although it has a more complex formula. Similarly, Z.  comprehensively demonstrated the balance between the model simulation accuracy and generalization ability, and then proposed a hybrid model that combines multiple harmonics with a linear function of LST-related factors.
As shown in Figure 26, based on the ATC model, DMLST can be used to estimate AMLST. Similarly, the ATC model can be extended to make it suitable for MMLST. Correspondingly, inputs and the annual period ω in Equations 17 and 18 become MMLST and 12, respectively. In addition, the idea of ABD can be further upgraded to average by month (ABM), namely, average by MMLST. The scheme first calculates MMLST using the methods mentioned in Section 5.4.5.2 and then retrieves AMLST by averaging all MMLSTs. Figure 27 illustrates the distribution of global AMLSTs in 2018 using the ABM-based method and MOD11A1/MYD11A1 products. Compared with Figure 24, upscaling from MMLST to AMLST, the data consistency of LST has been further improved. The ABM-based method fulfills the definition of AMLST well and is flexible because it does not need to fit the ATC model. However, the other two ATC-based methods can provide more parameters than ABM-based method except for AMLST.

Applications of LST
LST can provide unique information on thermal response to the variation in ET and SM, making it beneficial for the retrieval of such parameters. Additionally, LST has been frequently used to monitor agricultural drought, thermal environment, thermal anomaly, and climate change.

ET Estimation
LST provides an essential source of information that is closely linked to surface water status under the water-limited regime and therefore to the ET because varying SM conditions yield a distinctive thermal signature: moisture deficiencies in the root zone lead to vegetation stress and increase canopy temperatures, and depletion of water from the soil surface layer causes the soil to heat rapidly (Anderson et al., 2007;J. M. Chen & Liu, 2020;Pascolini-Campbell et al., 2021). Under the energy-limited regime where soil moisture is sufficient, LST and ET are both heavily dominated by atmospheric and radiative forcings. LST has been used for producing local,  Chen & Liu, 2020; R. L. Tang, Li, & Chen, 2011;Xu, Guo, et al., 2019;Xu, He, et al., 2019).

Absolute LST-Based Method
Using the absolute value of LST to replace or calculate aerodynamic temperature to estimate sensible heat flux with the bulk transfer equation in the single-source energy balance models is one of the most frequently used methods to estimate ET. In the absolute LST-based method, ET is usually obtained as the residual of SEB, and an excess resistance is generally incorporated, for example, in the well-known surface energy balance system (SEBS, X. L. Chen et al., 2021;Z. Su, 2002;Webster et al., 2017), the surface energy balance algorithm for land (SEBAL, Bastiaanssen, 2000;Bastiaanssen et al., 1998;Bhattarai et al., 2017;Liou et al., 2018;R. L. Tang, Li, Chen, et al., 2013), and the mapping ET model with internalized calibration (Allen, Tasumi, Morse, et al., 2007;Losgedaragh & Rahimzadegan., 2018;Tasumi, 2019).

Spatial Variability of LST-Based Method
Another approach when attempting to estimate ET under the condition of stressed soil moisture is using the spatial contextual information of LST to represent the SM availability or surface evaporative fraction in the simplified surface energy balance index (S-SEBI) model (Acharya & Sharma, 2021;Rahman & Zhang, 2019;Roerink et al., 2000), the Operational Simplified SEB model (Senay, 2018;Senay et al., 2014), and the LST versus vegetation index triangular or trapezoidal space ( Figure 28, Carlson, 2007;L. Jiang & Islam, 1999;Long & Singh, 2012;Minacapilli et al., 2016;R. L. Tang et al., 2010;R. L. Tang & Li, 2017a;R. Zhang et al., 2008). For the spatial variability of the LST-based method, ET (or EF) is estimated by interpolating the end-member values of ET (or EF); thus, much of the bias errors in estimating the aerodynamic temperature can be eliminated.

Viewing Angle Variability of LST-Based Method
The last method is to employ multi-angular measurements to retrieve soil/vegetation component temperatures from a mixed pixel LST and then separate soil evaporation/vegetation transpiration from ET within the framework of two-source energy balance models (Anderson et al., 1997;, 1999Song et al., 2020). These models compute surface fluxes without the measurement of air temperature , reduce the sensitivity of transpiration to input parameters, and improve the estimated ET (Song et al., 2020). Nevertheless, the subsequent dependence of ET estimation on the decomposed component temperatures increases (Chehbouni et al., 2001;. Regrettably, few studies have performed ET estimation by using the viewing angle variability of LST. Table 5 briefly summarizes the representative models for ET estimation using LST data (hereafter referred to as the LST-based ET models) and their advantages and disadvantages. Compared with the conductance-based ET models (e.g., Penman-Monteith model) that are independent of LST, the LST-based ET models better characterize the heterogeneity of SM at a fine spatial and temporal scale. Figure 29 displays three typical publicly available LST-based daily ET products generated by using the 5 km MODIS data, the 30 m Landsat 8 data, and the 70 m ECOSTRESS data. A summary of studies on the validation of LST-based ET methods against ground-based flux measurements shows an average RMSE value of approximately 50 W/m 2 and relative errors of 15%-30% in ET estimation (Kalma et al., 2008). Sensitivity analysis has shown that a 1 K variation in LST in a two-source model can incur a maximum bias of 29 and 30 W/m 2 with the average of 17 and 19 W/m 2 in sensible and latent heat fluxes, respectively (R. L. Tang, Li, Jia, et al., 2011;X. Zhan et al., 1996) and variations of 41% and 71% in sensible flux heat estimates in two single-source energy balance models (Kustas et al., 1989;Troufleau et al., 1995) have been observed as well when an uncertainty of 10% in LST (in °C) was varied (X. Zhan et al., 1996), suggesting that the improved LST retrieval could significantly decrease the uncertainty and error of ET estimates from the energy balance models.
Because LST is highly responsive to latent heat loss, the LST-based ET methods reveal obvious advantages for capturing ET spatial and temporal variations under water-limited conditions. However, several problems remain unresolved and require further investigation. The first problem is that these methods generally only provide instantaneous clear-sky surface ET at the satellite overpass time due to the instantaneity of LST retrieval from the remote sensing data (Wan & Dozier, 1996). To fulfill the requirement of ET in practice, the instantaneous ET estimated from LST must be further temporally upscaled to daily, weekly, or monthly ET (R. L. . However, upscaling is usually difficult because ET changes significantly in a diurnal cycle (R. L. Tang & Li, 2017b, 2017c. The second problem is the effect of the duration and amount of cloud cover, which may deteriorate the spatial and temporal continuity of LST to retrieve ET (He et al., 2020). Last, the accuracy of • Requires ground-based measurements as inputs • Highly sensitive to LST Z. Su (2002), , and Bastiaanssen et al. (1998) Use of spatial variability of LST S-SEBI, SSEBop, LST-VI triangle/trapezoid space models • Relatively insensitive to LST • Atmospheric conditions are usually required to be reasonably constant L. Jiang and Islam (1999), Roerink et al. (2000), Senay et al. (2014), R. L. Tang and Li (2017a) Use of the view angle variability N95, ALEXI, DisALEXI • Separates soil evaporation from vegetation transpiration • Requiring many groundbased measurements as inputs Anderson et al. (1997), , and Norman et al. (2003)  ET estimation from these models is usually subject to the observational angular anisotropy, leading to an inconsistency or inequality in the estimates of sensible heat flux and further ET.
In summary, the LST-based ET models developed since the early 1970s have been progressively refined and substantial progress has been made (Bouchra et al., 2020;J. M. Chen & Liu, 2020;Z.-L. Li et al., 2009). For the estimation of continuous ET with LST data, combining LST-based and conductance-based models may be useful through the integration of various remote sensing data because conductance-based models primarily use vegetation structural information derived from visible/near-infrared measurements (Garrigues et al., 2008;C. Jiang et al., 2017). For the angular anisotropy effect of LST, methodologies should be developed to improve the parameterization of ET models or accurately estimate the component temperatures of vegetation and soil from multispectral and multi-angular satellite measurements.

SM Estimation
Due to the effects of SM on the surface heating or cooling process under the soil, LST or LST variation is recognized to be sensitive to SM (Ghahremanloo et al., 2019). As a consequence, LST-based SM retrieval has become one of the most widely investigated research interests, with the increasing maturity of satellite LST retrieval Figure 29. Three typical publicly available LST-based ET products. (a) the MODIS daily ET product at 5 km spatial resolution (30 June 2017; https://data.tpdc.ac.cn/ zh-hans/data/df4005fb-9449-4760-8e8a-09727df9fe36/?q=EB-ET), estimated by the SEBS method; (b) the Landsat daily ET product at 30 m spatial resolution (15 July 2021; https://www.usgs.gov/core-science-systems/nli/landsat/landsat-provisional-actual-evapotranspiration), derived by the Operational Simplified Surface Energy Balance model (Senay, 2018;Senay et al., 2014), currently available only for the CONUS region; (c) the ECO3ETALEXI ET product at 70 m spatial resolution (24 May 2021; https://lpdaac.usgs.gov/products/eco3etalexiv001/) available for CONUS, derived by a physics-based SEB algorithm and the Atmosphere Land Exchange Inverse (ALEXI) Disaggregation algorithm (DisALEXI).
technology since the 1970s. Among all existing methods, SM retrieval based on TIR and VIS/NIR synergy has received more attention in recent decades (Z.-L. .

ATI-Based Method
Apparent TI (ATI) is one of the most widely used methods for estimating SM from the TIR and VIS/NIR synergy, which can be written as follows (Mitra & Majumdar, 2004): where ΔT is diurnal LST amplitude, α is the land surface albedo derived from VIS/NIR channels, and C is the solar correction factor that can be determined from latitude and solar declination for a given pixel.
The basic principle of the ATI method is that significant soil thermal conductivity contrast exists between wet and dry soils. It should be noticed that Equation 19 is only a simplified version of the ATI originally defined by Price (1985) where more surface physical parameters and meteorological factors were required. The primary motivation for the simplification of the ATI is to make it easier in practice with remote sensing data, because many TIR and VIS/NIR satellites provide ΔT and α with considerable accuracy. Figure 30 depicts the satellite-derived ATI and the synchronous Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) SM product over the southern part of the African continent on 26 March 2009. It notes that both ATI and SM have similar spatial distribution over the study region, indicating that ATI has significant potential to infer SM status. However, because the relationships between ATI and SM largely depend on soil texture, further investigation is required to transform ATI to volumetric SM content (S. Lu et al., 2009;Y. Lu et al., 2018;Nearing et al., 2012). Currently, two approaches have been widely used to obtain volumetric SM content from satellite-derived ATI: develop empirical relationships between SM and ATI, where a large amount of in situ SM measurements are commonly required (Paruta et al., 2021;Yuan et al., 2020), and use time series ATI to normalize ATI in a given study period under the assumption that the maximum and minimum ATI over a long time series can potentially refer to maximum (e.g., field capacity) and minimum (e.g., wilting point) SM content, respectively (Sohrabinia et al., 2014;Verstraeten et al., 2006).

Feature Space-Based Method
The triangular or trapezoid feature space that was constituted by the LST and fractional vegetation cover (FVC) in the two-dimension spatial distribution (denoted as LST-FVC) is another widely used method for estimating SM with various satellite data (Carlson et al., 1981;Goward et al., 1985;Moran et al., 1994). Explicitly, the basic hypothesis for the feature space method is that SM varies linearly with LST under a given vegetation condition, and LST is more sensitive to SM at low vegetation conditions than that at high vegetation coverage (Petropoulos et al., 2009). Figure 31 depicts a diagram of the trapezoid feature space. The four pixels (I, II, III, and IV) refer to dry bare soil, full vegetation cover with zero water availability, well-watered full vegetation cover, and saturated bare soil, respectively. Specifically, an additional hypothesis of comparing a trapezoid to triangular feature space is that LST continues increasing when FVC is at its maximum value under water stress conditions. For the feature space, SM availability (M 0 ) over the dry edge is considered to be zero, indicating that no water is available for evaporation or transpiration. In contrast, M 0 over the wet edge is equal to 1, showing a sufficient amount of water available for ET. As described, because SM is assumed to vary linearly with LST for a given FVC condition in the feature space, M 0 for any target pixel within the feature space can be written as follows: where A and B are the distances of the target pixel to the dry edge and wet edge for the same FVC, respectively, which can be determined by LST over the dry/wet edge and the target pixel. It notes that many investigations also determine the temperature vegetation dryness index (TVDI) from the feature space as TVDI = B/(A + B); thus, M 0 and TVDI have a mathematical relation that their sum is equal to 1 (Leng et al., 2019;Sandholt et al., 2002). The main distinction is that TVDI is only used in a remote sensing community, whereas M 0 is also a concept that frequently occurred in soil or agricultural sciences.
It notes that the spatial distribution of LST-and FVC-based-feature space generally requires homogeneous atmospheric conditions and varying vegetation coverage/SM status over the pixel window. Thus, it is probably subjective to determine the pixel window and dry/wet edges in practice. To this end, Leng et al. (2017) developed a novel concept of "pixel-to-pixel" feature space. Specifically, the pixel-to-pixel scheme of the feature space assumes a unique, virtual feature space exists corresponding to each satellite pixel over the given study area, and the theoretical dry/wet edge for any pixel can be exclusively determined by the atmospheric conditions over the pixel. Following the pixel-to-pixel scheme of the LST-FVC feature space, Leng et al. (2019) further developed a method to obtain all-weather SM in China (Figure 32). In this method, SM over MODIS clear-sky pixels was estimated with the pixel-to-pixel feature space, whereas SM over MODIS cloudy pixels was obtained via the disaggregation of coarse spatial resolution microwave-based SM products.

DTC-Based Method
Except for the aforementioned two widely used LST-based methods, two notable achievements of SM retrieval in recent years are the methods proposed by W. Zhao and Li (2013) and Leng et al. (2014) that take the advantage of the mid-morning period of the DTC information and the entire daytime period of the DTC information to estimate SM, respectively. The basic principle of the DTC-based SM retrieval methods by W. Zhao and Li (2013) and Leng et al. (2014) is that the DTC characteristic parameters represent the LST variation leading by energy per  unit, which largely depends on SM status. Figure 33 depicts the DTC curve and the linear relationship between LST and net surface shortwave radiation (NSSR) over the mid-morning period used in the method of W. Zhao and Li (2013), and Figure 34 depicts the elliptical relationship between DTC and the diurnal NSSR cycle used in Leng et al. (2014) method. It notes that for the elliptical relationship, the LST and NSSR data have been pre-processed to make them dimensionless. A major advantage of the DTC-based SM retrieval methods is that they obtain volumetric SM content directly; thus, the knowledge of soil texture and the process for decoupling volumetric SM content from satellite-derived SM proxies are not required. Based on the characteristic parameters in the DTC model and the elliptical relationship, W. Zhao and Li (2013)'s method and Leng et al.'s (2014) method can be written as follows: = 1 × 0 + 2 × 0 + 3 × + 4 × + 0 where T NM and t m are the LST rising rate normalized by the difference in the NSSR during the mid-morning and the time of the maximum LST occurs, respectively; x 0 , y 0 , a, and θ are the ellipse parameters derived from the elliptical relationship, namely, the ellipse center horizontal coordinate, ellipse center vertical coordinate, semi-major axis and rotation angle, respectively; and m i (i = 0, 1, 2) and n i (i = 0, 1, 2, 3, 4) are the model coefficients (m 3 /m 3 ) that can be obtained by LSM simulation.
Studies have attempted to estimate SM from the combination of TIR and microwave observations, which indicates a new direction for SM retrieval based on LST data (Amazirh et al., 2018;Mishra et al., 2018). Additionally, literature that investigates using LST and other auxiliary data to obtain high resolution SM from the microwave-based coarse SM products through disaggregation processes is increasing (J. Peng et al., 2017).
Microwave-based SM products are available, but no LST-based SM product is available. Hence, LST-based SM retrievals are generally validated in the stage of methodology developments. Most of the aforementioned LST-based SM retrieval methods exhibit an accuracy with RMSE of ∼0.06 m 3 /m 3 (Van Doninck et al., 2011;Verstraeten et al., 2006). In summary, most of the currently available LST-based SM retrieval has several major limitations including the requirement of decoupling volumetric SM content from LST-based SM proxies, the decreased retrieval accuracy over the dense vegetated conditions, and the undefined soil depth representativeness. Specifically, it should be noted that the ATI-based method and feature space-based method can only obtain a SM proxy (ATI or M 0 ) directly from satellite observations. As a consequence, satellite-derived SM proxies should be converted to volumetric SM content that desired in practices. In general, two approaches are commonly used to deal with the conversion: the one is to simply establish statistical relationships between satellite-derived SM proxies and volumetric SM content measurements; the other is through the assistance of soil texture. Compared to the ATI-based method and feature space-based method, the major advantage of the DTC-based method is the feasibility of estimating volumetric SM content directly without calibrations based on volumetric SM content measurements or soil texture. However, it requires meteorological data as auxiliary input. Reads are encouraged to obtain more information for these methods from a recent review article by Z.-L. . Besides, although a few literatures have highlighted the possibility of inferring deeper SM status using the aforementioned methods, it notes that the mainstream still considers satellite-derived SM as those contained within the top soil layer (e.g., ∼5 cm) (Babaeian et al., 2019;. Table 6 summaries the aforementioned LST-based SM retrieval methods in detail.

Agricultural Drought Monitoring
Due to the sensitivity of LST to SM and SM-associated agricultural drought conditions, LST-based indices have been developed to monitor agricultural drought worldwide. According to the data span used to form the indices, these agricultural drought indices can be generally classified into two major categories: (a) instantaneous LST-based indices; (b) time series LST-based indices.

Instantaneous LST-Based Indices
SM is commonly regarded as a measure of agricultural drought because it affects plant growth and productivity (Krishnan et al., 2006;Veihmeyer & Hendrickson, 1950). Hence, the aforementioned LST-based SM proxies such as the ATI and M 0 (TVDI) can be used as potential indicators for monitoring agricultural drought. Other widely used LST-based SM proxies that can be used as agricultural drought indicators include the vegetation supply water index (L. Li et al., 1998), crop water stress index (Jackson et al., 1981) and water deficit index (Moran et al., 1994). Table 7 shows a summary of these instantaneous LST-based agricultural drought indices and their advantages and disadvantages in detail.
Notably, the "instantaneous" used here is a relative term as the "time series" in the upcoming subsection, because many scholars have also used the 8 days composited MODIS LST to constitute the indices in practice. Another notable problem of such indices is that they are actually not commonly used to monitor agricultural drought but to infer SM conditions. This phenomenon is generally attributed to the following: these indices are determined from instantaneous or composited LST rather than using anomaly LST values over a long period. Hence, these indices are regarded as being spatially comparable to infer the relatively dry or wet conditions of soil within a specific region over the same scene of an image; however, they are not temporally comparable to describe the onset, duration, intensity, and end of the entire process of agricultural drought.

Time Series LST-Based Indices
Time series LST-based indices have been used more often to monitor agricultural droughts in the past few decades than instantaneous LST-based indices. The principle of the agricultural drought monitoring using time series LST data is that these LST-based indices were determined by calculating the extent of an anomaly or departure from the longer-term environmental baseline . The temperature condition index (TCI) is a simple time series LST-based index to determine temperature-related vegetation stress developed by Kogan (1995). The TCI is defined as where T s,max and T s,min are the maximum and minimum LST over a long time series, T s is the LST over the study period, and 100 is an enhance factor.
After the TCI, Kogan (2001Kogan ( , 2002 proposed a new index, namely, the Vegetation Health Index (VHI), to monitor agricultural drought from the combination of TCI and the Vegetation Condition Index (VCI). Specifically, the VHI integrates time series of LST and NDVI, which probably has an enhanced skill to monitor the stress of vegetation related to both temperature and water. The VHI is defined as where Θ and 1 − Θ refer to the relative contribution of the VCI and TCI, and NDVI max and NDVI min are the maximum and minimum NDVI over a long time series as the LST. In general, if moisture and temperature contribution during a vegetation cycle is unknown, Θ can be set as 0.5 in practice. In general, the VHI varies from 0 (extreme drought) to 100 (no drought), representing different dryness or wetness levels.
Except for the TCI and VHI that monitor droughts from the perspectives of temperature and water, a unique feature of ET for drought monitoring is that it describes both water availability and the rate at which it is consumed. Following the geostationary satellites-derived LST, Anderson et al. (2011) developed the Evaporative Stress Index (ESI) to measure anomalies in the ratio of actual to potential ET (PET) to monitor drought conditions. Specifically, the ratio of actual to potential ET can be written as where f PET is the ratio of actual to potential ET.
In the method proposed by Anderson et al. (2011), the standardized anomalies of f PET over a long period is expressed as a pseudo z score with a mean of 0 and a standard deviation of 1. The standardized anomalies were computed as PET is the averaged f PET over a long time series, and σ is the standard deviation.
The VHI and ESI remain among the most widely used satellite-derived indices for monitoring agricultural droughts by the US Drought Monitor to provide weekly drought assessments across the United States (https:// droughtmonitor.unl.edu/). Figures 35 and 36 depict the VHI and ESI on 16 November 2021, in the United States by referring to the report of weekly drought indices. Most of the western to central US regions probably suffer from various droughts from both the VHI and ESI during the investigated period.
In addition to these indices, several studies have used times series LST, NDVI, and other variables such as satellite-derived precipitation/SM products to monitor agricultural droughts. A typical composite index of this type is the Scaled Drought Condition Index (SDCI); it combines MODIS-derived LST and NDVI and the Tropical Rainfall Measuring Mission satellite-derived precipitation (Rhee et al., 2010). For the three variables in the determination of SDCI, each variable was scaled from 0 to 1 using the maximum and minimum values over a long time series to discriminate the effect of drought from normal conditions and then combined with the selected weights. Subsequently, a drought vector method was proposed from the anomaly indices by LST, NDVI, and the satellite-derived SM index (SWI), which was further used to quantify and predict drought in Northwest Africa during a 10-year period (Le Page & Zribi, 2019). In summary, unlike the instantaneous LST-based indices, the time series LST-based indices are temporally comparable, which is feasible for the quantification and monitoring of the beginning, development, and end of droughts for each pixel. However, comparing drought conditions for different geographic locations with the time series LST-based indices is difficult because the determination of the drought indices for a given pixel is independent of other pixels.
Due to the complicated process of drought events and various definitions of agricultural droughts, a critical issue with respect to agricultural drought monitoring is that most of the currently available methods still largely depend on mathematical statistics. Although many LST-based methods have been developed to monitor agricultural droughts at various spatiotemporal scales at present, most of these methods still lack of sufficient physical mechanisms, thus cannot well capture the accurate onset, development and end of real drought events.

Surface Urban Heat Island Monitoring
The phenomenon that urban LST is higher than that in suburban or rural regions is defined as the surface urban heat island (SUHI), which acts as a vital indicator to monitor the urban thermal environment (Oke, 1982). With different quantification methods and various LST sources, the global SUHI dynamics have been thoroughly investigated at multiple spatial and temporal scales (Arnfield, 2003;Oke, 1988;Rao, 1972;L. Zhao et al., 2014).
Generally, SUHI intensity (SUHII) is intuitively quantified by the LST difference between urban and rural regions as follows: In which the LST u and LST r denote the LST in urban and rural regions, respectively. It is prevalent in large-scale or long-term studies because of the computing convenience. However, the result is sensible to the identification of urban and rural regions (Clinton & Gong, 2013). The MODIS LST is widely used in revealing the spatial and temporal patterns of global SUHI because of its wide spatial coverage and long-term sequence (Chakraborty & Lee, 2019;S. Peng et al., 2012;Si et al., 2022;D. Zhou et al., 2019). Regarding to the SUHI at regional or national scales, the LST products with finer spatial resolution such as those derived from ASTER and Landsat increased in popularity, and the research has been showing a trend of exponential growth (D. . A global long-term SUHI was depicted, and a subtle increase trend of daytime SUHI at 0.03°C/decade was observed (Chakraborty & Lee, 2019). However, most of the current SUHI studies have used the simplified urban extent, which overlooked the urban expansion and may result in bias in SUHI quantification. Afterward, a dynamic urban extend was proposed to improve the depiction of the global long-term trend of SUHI, and Figure   At a city or local scale, the LST pattern-driven method can depict the LST distribution fields in certain urban clusters through complex spatial distribution models (e.g., Gaussian volume model) (Streutker, 2003). The magnitude of the regional LST and SUHI footprint can be detected to display the magnitude center of SUHI intensity and its spatial patterns (e.g., extension, inclination) simultaneously (Anniballe et al., 2014;Rajasekar & Weng, 2009). The SUHI footprint provides a superior approach to tracing the spatial pattern of SUHI, especially for local-scale studies. By using the Gaussian model, SUHI intensity generally shows an exponential decaying trend from the urban center to the rural area, and the footprint is beyond the boundary of the urban area in several studies (Q. Yang et al., 2019;D. Zhou et al., 2015). It is suitable for better monitoring of the SUHI, whereas the complex model fitting is inapplicable for global studies (Z.-L. . Notably, the spatial resolution of LST is essential to guarantee the accuracy of the SUHI footprint model. On this occasion, the LST product with finer resolution (e.g., Landsat, ASTER) is superior for the quantification of the SUHI footprint.
With respect to the temporal dynamic of SUHI, its long-term variations appear to be another essential issue. The linear regression was frequently used to characterize the interannual trends of SUHI, and several results have shown that the changing rate of SUHI significantly relies on surface biophysical characteristics, urban sprawl, and human activities (Yao, Wang, Huang, Liu, et al., 2021;W. Zhou et al., 2014). Additionally, non-parametric methods such as the Mann-Kendall test and Sen's slope have been employed to detect the existence of long-term trends and its changing rate (Mathew et al., 2016;Yao, Wang, Huang, Guo, et al., 2017). Recently, the ensemble empirical mode decomposition method derived from the digital signals field was introduced to determine the long-term variation patterns of LST.
As one of the most severe urban thermal environmental problems, abundant SUHI monitoring by LST has been increasingly developed. However, several issues should be especially concerned in the application. First, the precise quantification of SUHI intensity largely depends on the accurate identification of urban and rural regions, and SUHI studies have generally ignored the time-varing urban morphology accompanied by the urban sprawl.
To this end, the urban target should not be restricted to a fixed urban boundary, and the dynamic urban and rural regions should be accurately identified and updated, especially over rapidly-expanding regions. Second, the rural background with mixed land cover type exhibits heterogeneous LST components, which would result in incomparability among regions. Thus, investigating the effect of diverse rural backgrounds on SUHI research and improving the methodology of SUHI quantification for the comparison studies of cities are necessary (Si et al., 2022). Finally, concerning the spatial and temporal discontinuities of LST datasets, the frequency of valid LST observations will also disturb the integrality of SUHI for specific regions (e.g., global, regional, city scales) or over certain periods (e.g., annual, seasonal, monthly scales) as follows: (a) The spatial coverage of LST may be inconsistent and lead to SUHI being incomparable among regions. In a city with a discrete LST field, the reliability of SUHI quantification is reduced and should be specifically concerned; otherwise, it may further disturb the investigation of the SUHI pattern at larger scales; (b) the temporal discontinuities of LST datasets may impede long-term SUHI studies. In urban circumstances, the probability of cloudy days are generally higher due to the aggregation of urban buildings and facilities, and the insufficient frequency of valid LST observations disturbs the continuity of SUHI over the long term (Morris et al., 2001). Capturing long-term variations requires filling in the missing result either by the temporal upscaling of the LST data source or the temporal expanding of SUHI.

Thermal Anomaly Monitoring
When the LST changes for some reason and deviates from its normal distribution in spatial or temporal distribution, it usually called as LST Anomaly, and can be detected based on this principle from the remote sensing TIR data. Specifically, the positive (negative) thermal anomaly indicates that the observed LST was higher (lower) than the reference value. Normally, the LST values follow the order of Gaussian distribution or T-distribution in either the spatial or temporal dimension. Hence, the crucial issue in thermal anomaly monitoring is to determine a threshold to quantify the difference between the anomalous and normal pixels. A general function to define thermal anomaly with LST (LST_TA) can be written as follows: where ( ) refers to the LST of pixel (x, y) at time t, and and refer to the averaged LST and the standard deviation, respectively. The parameter a is a predetermined threshold to qualify the confidence interval of LST values, and usually can be set as greater than 2 with a confidence level over 95%. The pixel (x, y, t) can be regarded as an LST anomaly when LST_TA is greater than 0. In general, can be calculated from either the spatial or temporal or spatiotemporal dimension. Hence, thermal anomaly monitoring approaches can be divided into the spatial LST-based method, temporal LST-based method, and spatiotemporal LST-based method.

Spatial LST-Based Method
When the energy from the inner Earth reaches the land surface from the cracks, it often forms isolated and stable positive thermal anomalies, such as hot springs or coalfield fires. These anomalies do not change much over time and can thus be detected by the spatial LST-based method. For this situation, indicates the spatial averaged LST calculated from neighbor pixels within a sliding window at time t, which can be written as follows: where ( , , ) represents the LST of pixel (x i , y j ) at time t, the parameters m and n represent the sliding window size (m × n), and r refers to the number of valid LST pixels within the sliding window. Notably, anomaly events should have a certain spread (e.g., several gathered pixels), or they are most likely to be regarded as outliers with the spatial LST-based method. Specifically, the rule can be defined as Δ(LST_TA)>S, in which Δ(LST_TA) is the spread area of the spatially aggregated LST_TA pixels, and S is the predefined minimum threshold of spread area.
Geotherm and coalfield fire detections are two typical examples of thermal anomaly monitoring using the spatial LST-based method. Geotherm refers to the internal energy leaks out at the edge of tectonic plates and causes the LST to increase significantly (e.g., volcanic geotherms and hot springs). In general, a geotherm has a shallow depth and can increase LST; thus, it can be directly detected (Vaughan et al., 2012). In practice, there are two commonly used methods to detect a geotherm: using the D/N LST difference instead of daytime LST or nighttime LST in Equation 29 (i.e., day-night comparison method) or taking advantage of LST and many other ancillary data (comprehensive analysis method) (K. Wang et al., 2019). In practice, the comprehensive analysis method has received more attention. A typical example is the extraction of volcanic geotherm over the Main Ethiopian Rift region by Yosef et al. (2019). Normally, the temperature would consistently increase approximately 3 K for each 100 m of depth. As in the volcanic geotherm region, much more heat flux and significantly higher LST can be observed by satellite sensors than other regions (DiPippo, 2005;Q. Qin et al., 2011). Based on bands 10 and 11 of the Landsat 8 satellite, Yosef et al. (2019) calculated the LST by using an SW algorithm; next, pixels with a higher LST (3-9 K) than background value were primarily selected as potential heat sources (Figure 38a) based the spatial windows and the aforementioned rule of Q. Qin et al. (2011). Subsequently, by combining the LST-derived potential heat sources and ancillary data of a geological map, a topographic map, and DEM, the volcanic geotherm was finally determined (Figure 38b).
Coal is a fuel that is subject to catching fire by spontaneous combustion, a nearby forest fire, and improper mining technology. Because LST is higher in the zones of underground coal fire than in their surrounding areas, there is a possibility to monitor the coal fire by using the spatial LST-based method (Gangopadhyay et al., 2006).
Conventionally, the spatial LST-based method has been widely used in coal fire detection; however, it merely approximately indicates the segment threshold of temperature for coal fires. More specifically, the threshold can be set with the confidence level in accordance with the field survey of the LST anomaly in typical coal fire regions (Gangopadhyay et al., 2006;Xue et al., 2018). In addition to the spatial pattern and size, the spreading direction and development of dynamic coal fire can be captured by employing multi-temporal satellite LST, which is essential for controlling coal fire propagation and preventing energy loss (Biswal & Gorai, 2020;Huo et al., 2015). Notably, nighttime LST presents more accurate results than daytime LST for coal fire monitoring due to the exclusion of solar radiation at nighttime (Biswal et al., 2019). Figure 39 presents the coal fire propagation from 2006 to 2015 in the Jharia coalfield in India (Singh et al., 2021).

Temporal LST-Based Method
Thermal anomaly events are caused by the dramatic changes in meteorology that may change their spread over time, and due to the heat exchange between land surface and air, the change of the skin LST will be consistent with the air temperature. In this situation, the thermal anomaly can be detected by the temporal LST-based method, and indicates the temporal averaged LST calculated from LST time series for pixel (x, y), which can be written as follows: where ( ) represents the LST of pixel (x, y) at time t k , p represents the LST time series, and r represents the number of the valid pixels. Notably, the anomaly events should last for a range of time (e.g., several days) in a certain period; otherwise, they are most likely to be regarded as outliers when using the temporal LST-based method. The rule can be defined as Ω(LST_TA)>TT, where Ω(LST_TA) is the last time of LST_TA for the same pixel, and TT is the predefined minimum threshold of last time.
Studies have indicated that the rise of global LST has increased (decreased) the frequency of positive (negative) anomaly events in climate. In general, a heat wave is a classic positive thermal anomaly, which has obvious temporal characteristics including the onset, duration, intensity, and end. Figure 40 shows the heat wave of Eastern Europe in 2010 detected using the temporal-based method with the MxD11A1 LST product from 2000 to 2020. The statistical LST-derived indices (e.g., heatwave number, longest duration, total days, average magnitude, and amplitude-average) have also been used to evaluate heatwave activities (Perkins & Alexander, 2013). Several physical models have also been proposed to investigate the causes, developments, and impacts of heatwaves by using satellite-derived LST products (Schumacher et al., 2019). In contrast with a heatwave, a negative thermal anomaly usually refers to an unexpected rapid decrease in LST over a short time (i.e., 24 hr or less). In general, nighttime LST is frequently used to monitor a negative anomaly, because of the weak impact of wind on LST under such conditions.

Spatiotemporal LST-Based Method
When the surface environment (e.g., earthquake zone or volcano) changes, the stable LST may also significantly changed in both spatial and temporal dimensions. For such thermal anomaly detection, spatiotemporal LST is commonly required. Thus, the can be written as follows: where LST (x i , y j , t k ) represents the LST of pixel (x i , y j ) at time t k , and r represents the number of the valid num of the pixels within the sliding window in the time series.
Volcano is a procedure in which magma of the upper mantle rises to the surface along cracks in the crust and forms mounds of geological bodies, which can significantly change the temperature of land surface or air. The use of infrared data to detect volcano can be seen since 1960s (Wooster & Rothery, 2000), and remote sensing infrared data monitoring can be seen in 2004 (Pieri & Abrams, 2004). Due to that the volcanism is much more strenuous than the other thermal anomaly activates, many simple methods have been used to extract it, such as threshold-based method or context-based method. In fact, volcanoes are stable for a long time, and they will produce a lot of soot or water vapor when they become active. The recent studies have pointed that the volcanism will lead to intense surface warming accompanied by many volcanic soot or water vapor (Joshi & Jones, 2009;McCormick et al., 1995;Millán et al., 2022). Based on the above method, Pavlidou et al. (2017) successfully detected subtle spatiotemporal signal fluctuations to monitor volcanic activity of Mount Etna and Virunga National Park (Figure 41).
Thermal anomaly detection is an important application of LST products and has been developed for decades. Nevertheless, there remain many limitations to overcome. The primary issue is that most of the current applications rely on moderate-spatial/high-temporal LST products (e.g., MODIS and AVHRR). However, these LST products are not suitable for the detection of an anomaly with a relatively small size. In addition, almost all thermal anomaly studies are inevitably affected by clouds or cloud shadows; thus, how to remove their influence also requires further investigation. Last, most of the currently available methods for monitoring thermal anomaly primarily take the advantage of LST statistics, making the core issues of the anomaly detection (e.g., determination of the normal status and threshold) somewhat depend on the availability of LST products. As a consequence, the thermal anomaly detected by LST can potentially correspond to different anomalous issues or refer to a pretended anomaly event. Hence, physical models are encouraged to monitor thermal anomalies in further developments.

Climate Change Indices
In general, a small shift or variability in mean temperatures can lead to a disproportionally larger change in extreme temperatures (Caloiero, 2017;Field et al., 2012), indicating that extremes encapsulate more notable indicative signals for climate change. Thus, monitoring extremes is very important in the climate community. On the basis of the daily maximum and minimum temperatures, the Expert Team on Climate Change Detection and Indices developed 16 climate change indices for measuring extreme temperatures (L. V. Alexander et al., 2006). These widely used indices (Table 8) cover all aspects of extreme temperatures, including the intensity, range, frequency, and duration. These indices have mainly been developed based on air temperature from site observations or reanalysis data, which are either unevenly distributed globally or have coarse spatial resolution. Therefore, using satellite-derived instantaneous LST observations to estimate extreme indices improve climate change detection. However, due to the time discontinuity and limited historical data of current products, only intensity and range indices have been studied.

Intensity Index of LST
According to the definition, the precondition of retrieving intensity LST indices (TXx, TNx, TXn, and TNn) is to obtain the daily maximum and minimum LST. As aforementioned, the overpass times of afternoon satellites are close to the occurrence of these two extreme LSTs. Because the maximum LST usually occurs in clear-sky conditions, using instantaneous LST at 13:30 as the proxy of daily maximum LST is reasonable. When TXx and TXn are discussed, this rationale would strengthen because the temporal aggregation can remove the uncertainty from instantaneous retrieval to some degree. Correspondingly, Mildrexler et al. (2011) revealed the spatial patterns of the thermal pole of Earth by using MYD11C2 products, Azarderakhsh et al. (2020) analyzed the temporal variations in extreme LSTs across the hottest place on Earth based on MYD11A1 products, and Mildrexler et al. (2018) explored correlations between extreme LSTs from MYD11A2 and other events, such as land cover change and drought. Figure 42 illustrates an example that annual maximum LSTs are closely related to land cover types ( Figure 42a) and that their density probability changes effectively reflect the spatial variation in corresponding land cover (Figures 42b and 42c). However, saturation temperatures of satellite sensors (Prata, 2000) and active fire (Azarderakhsh et al., 2020)

Range Index of LST
Two methods have been used to estimate diurnal temperature range (DTR). One method is based on MDC (Aires et al., 2004;Holmes et al., 2015;Ignatov & Gutman, 1999;Norouzi et al., 2015;D. Sun et al., 2006) or ATC models (Y. Yu et al., 2021;W. Zhao, He, Yin, et al., 2019) in which the amplitude parameter represents DTR. The other method is relatively simple, namely, using day-night temperature differences of afternoon satellites to describe DTR (Ruzmaikin et al., 2017;Sharifnezhadazizi et al., 2019;Stuart-Menteth et al., 2003). The focus of these studies is on two topics. The first topic is exploring how to improve the reconstruction of temperature cycles, such as analyzing the requirement and distribution of measurements used (Ignatov & Gutman, 1999), developing new temporal interpolation algorithms for the situation when only a few LST measurements are available (Aires et al., 2004), and improving the adaptability of ATC models to simulate multi-year LST change (Xing et al., 2020). By contrast, the second topic emphasizes analyzing the spatial-temporal pattern of DTR (Ruzmaikin et al., 2017;Y. Yu et al., 2021) and its relationship with land surface conditions in response to climate change (Sharifnezhadazizi et al., 2019;D. Sun et al., 2006). Figure 43 illustrates the global distribution of monthly DTRs in 2018-01 and 2018-07 using MYD11A1 products. First, DTRs in high-latitude zones are usually larger than those in low-latitude zones. In addition, DTRs in the northern hemisphere are more significant in summer (i.e., 2018-07) than in other seasons, whereas DTRs in the southern hemisphere are more pronounced in winter (i.e., 2018-01) than in other seasons. Notably, the DTC model and ABD scheme are also applied to estimate DTC at a monthly scale. Specifically, daily DTR is first derived from the calibrated DTC model; next, all available daily DTRs are averaged to obtain monthly DTR. The ABM scheme can also be used to calculate annual DTR.

Perspectives
Accurately acquiring LSTs at different spatial and temporal scales is crucial for understanding the Earth's surface water and energy balances, such as evapotranspiration, drought, thermal anomalies, and climate change, in various studies. A comprehensive analysis and discussion of potential development directions and perspectives on LST retrieval and validation methods, products, and applications are vital for professional TIR communities and LST product users. As stated in the review of satellite-derived LST by Z.-L. Li et al. (2013), several approaches have been suggested for improving LST estimation from space-based measurements. In the last decade, some advances have been made in those directions, as summarized in Sections 3-5 of this paper. However, the problems associated with these directions have not been completely or effectively solved. Further investigations are required to reduce the gaps between the status of LST product development and the requirements of LST applications in various fields. Furthermore, new challenges regarding the generation and applications of satellite-derived LST products have been identified owing to the advances made in the last decade, which will require more attention in the near future.

Generating Instantaneous LST With Spatial Continuity and Temporal/Angular Consistency From Satellite Observations Under All-Weather Condition
TIR remote sensing is the best way to obtain LST with a high retrieval accuracy at regional and global scales. However, its LST applications in various fields is severely limited by its inability to penetrate clouds. Many methods (statistical regression, spatiotemporal reconstruction, and PMW-based methods) to obtain spatially continuous LST products have been proposed in recent decades. As discussed in Section 5.4.1, although progress has been made in reconstructing the LST of cloud-covered pixels, the proposed methods still suffer from low accuracy and poor practicability. The merging of TIR-and PMW-based LST is regarded as the most promising option to develop an all-weather LST product at the current technical level, considering the complementary properties of TIR-and PMW-derived LST. However, new efforts are still required to address the issues of their surface penetration depth discrepancy as well as the low accuracy and coarse resolution of PMW-based LST estimations.
In addition to the impact of cloud cover, LST varies strongly with both observation time and angle. In previous studies, as mentioned in Sections 5.4.2 and 5.4.3, most existing methods have been developed to separately correct the temporal or angular effects of LST products. The temporal and angular effects of LST products are coupled and should be corrected together. Therefore, a robust and effective temporal-and angular-normalized model should be developed for the concurrent correction of the two effects. The incorporation of different time and angle information acquired by multi-source LSTs from different satellite platforms could be a possible solution to drive this joint correction.

Developing Application-Oriented LST Characteristic Quantity Products
Mean and extreme LSTs find numerous applications in fields such as meteorology and climate change research (Table 7) where instantaneous LSTs do not, even though they have continuous spatial coverage and consistent observation time/angles. Therefore, it is urgently necessary to develop new application-oriented products for these LST characteristic quantities. As mentioned in Sections 5.4.5 and 6.6, some methods have been proposed to estimate the daily, monthly, and yearly mean LSTs and intensity/range indices of LST. However, their accuracies are modest (no better than 1.5 K). Furthermore, current methods are mainly empirical (e.g., TSA and MWA models) or derived from ideal conditions (e.g., DTC and MDC models), and therefore cannot be used for all applications. Therefore, developing an improved model with a satisfactory accuracy and clear physical meaning for estimating the temporal-scaling LST characteristic quantity is an important direction for future research. Until now, the frequency indices (i.e., TN10p, TN90p, TX10p, and TX90p) and duration indices (WSDI and CSDI) of LST have been lacking as the calculation of these percentile-based indices requires historical data (from 1961 to 1990) as the reference (L. V. Alexander et al., 2006). Compensating for this data gap is also a challenge. A possible solution is to develop the relationship between satellite-derived LST and other variables with the historical record (such as reanalysis data) and extend this relationship to the period without satellite LST. Alternatively, one can use existing data rather than historical data as the reference and utilize optimization methods, such as a bootstrap procedure (X. , to reduce the uncertainty caused by the limited record.

Constructing a Unified Archiving System for LST Products Storing and Distribution
From instantaneous LST products acquired from different satellite observations to application-oriented LST products for different scopes, the emergence of various LST products at different levels brings an urgent need to gather them across the globe. Therefore, a global or regional archive and data portal that facilitates data access would significantly promote progress in TIR remote sensing and its applications.
To refine the use of these LST products, a systematic standard should be established and applied by building a uniform level system with the definition of the product levels, establishing a uniform assessment criterion to label the quality of each product, and unifying the naming system of different LST products. Currently, although there is some progress in the definition of satellite-derived LST product levels, as discussed in Section 5.1, there remains a lack of consideration for application-oriented LST products. Therefore, the construction of product portals and standards will enable an improved understanding of LST products and their specific applications.

Establishing a Prototype for Validating LST Products
Long-term, temporal-, angular-consistent, and all-weather LST products would greatly promote LST applications in various fields on regional and global scales. Prior to these applications, robust and comprehensive validations are essential to evaluate the accuracy of LST products with different temporal and spatial resolutions. According to the review of current validation studies in Section 4, current LST product validation is mainly conducted with a small group/team over a homogenous surface. However, a global and persuasive LST product requires accurate evaluation of LST products with different temporal and spatial resolutions. Therefore, three suggestions are recommended for LST product validation in the near future: 1. Development of new LST validation methods for heterogeneous and non-isothermal surfaces. Combining ground-based multipoint measurements with high spatial-resolution LST acquired from unmanned aerial vehicle (UAV) platforms or low-orbit satellites could be a possible way to validate LST products over heterogeneous and non-isothermal surfaces using a level-by-level validation strategy. 2. Establishing a comprehensive validation scheme using multiple validation methods based on the availability of reference data and considering the representativeness of land cover type and vegetation phenology on a global scale. 3. Constructing a joint network of different research teams and institutions around the world to collect and release in situ LST measurements with a unified operating specification over well-maintained or newly installed sites. A referable example is the Land Product Validation (LPV) subgroup of the Committee Earth Observing Satellites' Working Group on Calibration and Validation (CEOS WGCV), which pursues the coordinated validation of land products using standardized protocols.

Analyzing Spatial-Temporal Patterns and Drivers of Global Long-Term LSTs
There is consensus that temperature variations are closely related to climate change. Currently, global warming is studied mainly with the help of surface air temperature measured 2 m above ground level. Compared to station-based air temperature, satellite-derived LST may better reflect the heat properties of the land surface (Mildrexler et al., 2011), and can provide uniform and complete spatial coverage with a higher spatial resolution. When the problems mentioned in Section 5.4 are resolved, it will be possible to explore the spatiotemporal patterns of long-term LST variations across the entire planet. Some preliminary studies have analyzed the linear trend of inter-annual LST changes within a relatively short period (e.g., less than 20 years) using traditional methods such as ordinary least squares (Susskind et al., 2019) and the Mann-Kendall test (Xing et al., 2020;Y. Yu et al., 2021). The gradual trend of multiyear LST may be piecewise linear or nonlinear. In addition, intra-annual variations (such as seasonal characteristics) and abrupt changes (such as those induced by fire disasters) can provide crucial indicator signals. Many methods have been proposed to capture the above three types of changes (i.e., inter-annual, intra-annual, and abrupt changes) in vegetation index or reflectance time series (Jamali et al., 2015;Verbesselt et al., 2010;K. Zhao et al., 2019); however, there is less concern regarding LST change detection. Therefore, it is necessary to improve and upgrade these methods so that they can be used to reveal LST variations at different temporal (such as intra-and inter-annual) and spatial scales (such as latitude bands, climate zones, administrative areas, and the globe).
Once LST variations have been characterized, the next task is to understand the mechanism of changes and identify their attributions, as this information is helpful for formulating climate mitigation policies. As LST is directly regulated by the radiative forcing and energy exchange efficiency between the land surface and the atmosphere, the derived LST variation could be attributed to the perturbation of other energy budget terms with the help of the SEB equation (Juang et al., 2007;Lee et al., 2011;. This method can be fully implemented based on observational data, explaining LST variation from the perspective of local biophysical effects. However, non-local biogeochemical effects play a predominant role in LST variation (Perugini et al., 2017). Accordingly, Earth system models (ESM) (C.  and advanced statistical models (Z. Y.  can explore all the driving forces of LST variations. It should be noted that large uncertainties exist in ESM simulations (Duveiller et al., 2018) and statistical models lack clear physical meaning. Therefore, a combination of model simulations and observation-driven methods may be a promising method to attribute LST variations.