Combining Spatiotemporally Global and Local Interpolations Improves Modeling of Annual Land Surface Temperature Cycles

: Annual temperature cycle (ATC) models are widely used to characterize temporally continuous land surface temperature (LST) dynamics within an annual cycle. However, the existing ATC models ignore the spatiotemporally local correlations among adjacent LST pixels and are inadequate for capturing the complex relationships between LSTs and LST-related descriptors. To address these issues, we propose an improved ATC model (termed the ATC_GL), which combines both the spatiotemporally global and local interpolations. Using the random forest (RF) algorithm, the ATC_GL model quantiﬁes the complex relationships between LSTs and LST-related descriptors such as the surface air temperature, normalized difference vegetation index, and digital elevation model. The performances of the ATC_GL and several extensively used LST reconstruction methods were compared under both clear-sky and overcast conditions. In the scenario with randomly missing LSTs, the accuracy of the ATC_GL was 2.3 K and 3.1 K higher than that of the ATCE (the enhanced ATC model) and the ATCO (the original ATC model), respectively. In the scenario with LST gaps of various sizes, the ATC_GL maintained the highest accuracy and was less sensitive to gap size when compared with the ATCH (the hybrid ATC model), Kriging interpolation, RSDAST (Remotely Sensed Daily Land Surface Temperature), and HIT (Hybrid Interpolation Technique). In the scenario of overcast conditions, the accuracy of the ATC_GL was 1.0 K higher than that of other LST reconstruction methods. The ATC_GL enriches the ATC model family and provides enhanced performance for generating spatiotemporally seamless LST products with high accuracy.


Introduction
Land surface temperature (LST) is one of the key parameters for determining landatmosphere energy exchanges and plays a crucial role in various applications, such as the monitoring of the urban thermal environment [1,2], estimation of soil moisture and evapotranspiration [3], and detection of forest fires [4,5]. Satellite thermal remote sensing provides a straightforward and consistent method for obtaining regional LSTs. However, due to cloud contamination and sampling error, more than half of the satellite-derived LST image pixels are invalid or entirely absent [6], and thus, the associated potential applications are severely limited.

Data
This study employed satellite data (MODIS LST, NDVI, and land cover type), in situ measurements (SAT), and auxiliary data (spatiotemporally seamless all-weather LSTs (i.e., the LSTMW data) and DEM) ( Table 2). Table 2. Features of the data and products used in this study.  * The LST missing rate is defined as the proportion of the missing LSTs within a yearly cycle in a certain region.

Data
This study employed satellite data (MODIS LST, NDVI, and land cover type), in situ measurements (SAT), and auxiliary data (spatiotemporally seamless all-weather LSTs (i.e., the LST MW data) and DEM) ( Table 2). The present study used the daily LST product (MOD11A1/MYD11A1, with a spatial resolution of 1 km) [44,45], 16-day NDVI product (MOD13A2, 1 km), and yearly landcover-type product (MCD12Q1, 0.5 km) ( Table 2). These products were obtained from NASA Earth Science Data. Daytime LSTs from MOD11A1 (~10: 30) in 2018 were used for LST reconstruction and model evaluation under clear-sky conditions (Strategy 1 & Strategy 2). The quality control (QC) of MODIS products was used to identify clear-sky or overcast conditions. The pixel was under cloud-free conditions when the corresponding QC was labelled as ' 10' or '11'. The NDVI product MOD13A2 (Section 3.1.2) and land cover type product MCD12Q1 (Section 4.1.2) in 2018 were used as model inputs and evaluations over different land cover types, respectively. The clear-sky LSTs of LST MW and the corresponding NDVI of MOD13A2 were used to reconstruct the LSTs under overcast conditions (Section 3.2). Note that the spatial and temporal resolutions of MCD12Q1 (0.5 km) and MOD13A2 (16-day) were resampled to 1 km and 1 day, respectively, to match those of the LST product by linear interpolation.

In Situ Measurements
In situ daily maximum and minimum SATs were derived from 2479 meteorological stations in China from 2018. To match the resolution of the LST data, the SATs were resampled to raster images with a resolution of 1 km after considering the impact of evaluation on temperature [22,46]. The daily mean SAT calculated as the average of daily maximum and minimum temperature was used to reflect the daily SAT fluctuations [46]. Subsequently, we calculated the daily mean SATs of each pixel [22,46] and used these values to quantify daily LST fluctuations.
In addition, in situ measurements obtained from seven Surface Radiation Budget Network (SURFRAD) sites in 2018 were used to evaluate the ATC_GL under overcast conditions ( Figure 2). We calculated the in situ LST from the upward and downward longwave radiation based on the Stefan-Boltzmann law [10,47], which can be expressed as follows: where R Long_up and R Long_down represent the upward and downward longwave radiation, respectively; σ denotes the constant of Stefan-Boltzmann (σ = 5.67 × 10 −8 W·m −2 ·K −4 ); ε b denotes broadband emissivity estimated from MODIS; ε 31 and ε 32 represent narrowband emissivity in MODIS Channels 31 and 32, respectively.

Auxiliary Data
The LSTMW data were used to evaluate the performance of the ATC_GL under overcast conditions. The LSTMW data were generated by merging thermal and microwave observations with the spatial and temporal resolutions of 1 km and 1 day, respectively [13]. The accuracy of LSTMW is about 1.5 K and there is no significant difference in accuracy under clear-sky and overcast conditions [13]. The LSTMW has two advantages over in situ LSTs: (1) the spatial and temporal resolutions of LSTMW are consistent with those of MODIS LST, and therefore, errors in LST retrieval and scale differences between satellitederived and ground-based LSTs are reduced; and (2) the available pixel-based LSTMW data are significantly larger than the ground-based LST data so that validations should be more representative [18,35,48]. However, the LSTMW only covers the northeast of China due to some limitations of microwave data (e.g., limited data coverage). Therefore, the LSTMW was used as the reference value for evaluating model performances under overcast conditions.
The yearly DEM product from SRTM3 was included in the quantification of daily LST fluctuations. The DEM product, with a spatial resolution of 90 m, was resampled to 1 km to match the resolution of the LST product.

Methodology
This section describes the framework of both model construction and evaluation (Figure 3). The construction of the ATC_GL involves (1) the temporally global reconstruction of annual LST dynamics (Section 3.1.1) and (2) spatiotemporally local estimation of daily LST fluctuations (Section 3.1.2). In terms of model evaluation, the performance of the

Auxiliary Data
The LST MW data were used to evaluate the performance of the ATC_GL under overcast conditions. The LST MW data were generated by merging thermal and microwave observations with the spatial and temporal resolutions of 1 km and 1 day, respectively [13]. The accuracy of LST MW is about 1.5 K and there is no significant difference in accuracy under clear-sky and overcast conditions [13]. The LST MW has two advantages over in situ LSTs: (1) the spatial and temporal resolutions of LST MW are consistent with those of MODIS LST, and therefore, errors in LST retrieval and scale differences between satellite-derived and ground-based LSTs are reduced; and (2) the available pixel-based LST MW data are significantly larger than the ground-based LST data so that validations should be more representative [18,35,48]. However, the LST MW only covers the northeast of China due to some limitations of microwave data (e.g., limited data coverage). Therefore, the LST MW was used as the reference value for evaluating model performances under overcast conditions. The yearly DEM product from SRTM3 was included in the quantification of daily LST fluctuations. The DEM product, with a spatial resolution of 90 m, was resampled to 1 km to match the resolution of the LST product.

Methodology
This section describes the framework of both model construction and evaluation ( Figure 3). The construction of the ATC_GL involves (1) the temporally global reconstruction of annual LST dynamics (Section 3.1.1) and (2) spatiotemporally local estimation of daily LST fluctuations (Section 3.1.2). In terms of model evaluation, the performance of the ATC_GL was compared with that of several popular LST gap-filling methods in three scenarios: (1) with randomly missing LSTs, (2) with LST gaps of various sizes under clear-sky conditions, and (3) under overcast conditions (refer to Section 3.2). In order to help readers better understand the mathematical symbols, we have listed their associated meaning in Table 3.

Abbreviations Descriptions Ts
The modeled LST on day d relative to the spring equinox Tg The general seasonal cycle of the intra-annual LST dynamics ΔTs The daily LST fluctuations gATCT The Tg modeled by the ATCT fRF The daily LST fluctuations modeled using the RF algorithm X The synoptic and surface parameters Ts -The annual mean LST

A1, A2
The two amplitudes of the ATCT ω1, ω2 The two constants calculated as 2πN -1 and 4πN -1 (N = 365) θ1, θ2 The corresponding phase shifts relative to the spring equinox gATCT_a The seasonal cycle in the intra-annual SAT dynamics determined by the ATCT Figure 3. The flowchart of the ATC_GL and the evaluation strategies in different scenarios. T g denotes the intra-annual LST cycle dynamics, ∆T s and ∆T a denote the daily fluctuations in LSTs and SATs, respectively, and LST, SAT, NDVI, and DEM denote land surface temperature, surface air temperature, normalized difference vegetation index, and digital elevation model, respectively. Table 3. All the mathematical symbols and their meaning in the following formulas.

Abbreviations Descriptions
T s The modeled LST on day d relative to the spring equinox T g The general seasonal cycle of the intra-annual LST dynamics ∆T s The daily LST fluctuations g ATCT The T g modeled by the ATCT f RF The daily LST fluctuations modeled using the RF algorithm X The synoptic and surface parameters T s The annual mean LST The two amplitudes of the ATCT ω 1 , ω 2 The two constants calculated as 2πN −1 and 4πN −1 (N = 365) θ 1 , θ 2 The corresponding phase shifts relative to the spring equinox g ATCT_a The seasonal cycle in the intra-annual SAT dynamics determined by the ATCT The two ATCT amplitudes of SAT

Modeling of Intra-Annual Land Surface Temperature Dynamics by Combining Spatiotemporally Global and Local Interpolations
Intra-annual LST dynamics can be mainly decomposed into two temporal categories: (1) the general seasonal (or monthly) cycle responding to solar radiation, and (2) the daily fluctuations caused by variations in synoptic conditions and surface thermal properties [19,22]. The ATC_GL addresses these categories separately. At the temporally global (seasonal) scale, the general seasonal LST dynamics (T g ) are modeled directly with the ATCT (see Section 3.1.1), a frequently used ATC model [17,19]. At the spatiotemporally local scale, the daily LST fluctuations (i.e., ∆T s ) are statistically modeled using the RF algorithm and a series of LST-related descriptors related to synoptic and surface property variations (see Section 3.1.2). The general formula for the ATC_GL is written as follows: where T s is the modeled LST on day d relative to the spring equinox; T g and ∆T s are the general seasonal cycles of the intra-annual LST dynamics and the daily LST fluctuations, respectively; T g and ∆T s are the intra-annual LST cycle dynamics and the daily LST fluctuations, respectively; g ATCT denotes the ATCT and is used to reconstruct the T g ; f RF represents the RF algorithm and is used to capture daily LST fluctuations through a series of synoptic and surface parameters (X), such as SAT, NDVI, DEM, and so on.

Reconstruction of Seasonal LST Cycle at the Temporally Global Scale
The existing ATC models are temporally global interpolations that capture the annual LST cycle by using sinusoidal functions to model periodical variations of solar irradiation [17,31]. Previous studies have demonstrated that the accuracy of ATCT modeled by a second harmonic function is generally better than the ATCO equipped with only a sine function, mostly because the annual cycle of solar irradiation is characterized by two annual maxima instead of one at low latitude regions [17,31]. To reconstruct seasonal LST cycles (T g ), we employed the ATCT, which can be expressed as follows: where g ATCT denotes the modeled LST on day d relative to the spring equinox, T s is the annual mean LST, A 1 and A 2 are the ATC amplitudes, θ 1 and θ 2 are the corresponding phase shifts relative to the spring equinox, and ω 1 and ω 2 are two constants calculated as 2πN −1 and 4πN −1 , respectively, where N is the number of days in an annual cycle (N = 365). The available ATC models can describe the variation of LST within an annual cycle but are not suitable for characterizing short-term LST fluctuations, e.g., on the daily scale. The differences between the observed and fitted LSTs (denoted as the daily LST fluctuations) can be expressed by Equation (4): In this study, we incorporated LST-related descriptors and the RF algorithm to simulate the daily LST fluctuations (∆T s ) at the spatiotemporally local scale. The number of decision trees and predictors of each split node were determined by the grid search method [48]. ∆T s is calculated using the following formula: The relationships between ∆T s and the LST-related descriptors are determined by three factors: (1) the type of LST-related descriptors, (2) the type of regression tool (or statistical algorithm), and (3) the spatiotemporal window used for regression. For the type of LSTrelated descriptors, ∆T s is closely related to surface parameters, such as NDVI and DEM, and is also synchronized with SAT under different synoptic conditions [11,19,22]; thus, we incorporated SAT, NDVI, and DEM to describe ∆T s . However, we did not incorporate SAT directly but included its daily variations (i.e., ∆T a ), which can be estimated with the following formula: where T a denotes the SAT, ∆T a denotes its daily fluctuations, and g ATCT_a denotes the seasonal cycle in the intra-annual SAT dynamics determined by the ATCT. A 3 and A 4 are the ATC amplitudes of SAT; θ 1 and θ 2 are the corresponding phase shifts relative to the spring equinox of SAT. For the regression tool, we used the RF algorithm to model the regression between ∆T s and the LST-related descriptors. The RF algorithm shows especially high performance in modeling complex relationships among parameters [49] and is widely used in remote sensing communities [15,50]. For the spatiotemporal window used for regression, we employed a spatiotemporally adaptive local window to incorporate valid (i.e., cloud-free) pixels. The initial window size for the spatial dimension was set to 9 × 9 km 2 and that for the temporal dimension was set to 1 day. The spatiotemporal window was expanded when the number of valid LSTs was lower than the threshold. Starting from the initial size, the spatiotemporal window size was iteratively increased by 2 pixels and 2 days from the spatial and temporal dimensions, respectively. Further discussions on the setting of the spatiotemporal window and the increased sizes of spatial and temporal windows are provided in Figures A1 and A2 and Table A1 of Appendix A. We further repeated this process until there were sufficient valid pixels included. The threshold of the valid pixel number was set to 10 [18]. Notably, the proposed ATC_GL predicts each missing pixel separately during the reconstruction process.

Evaluation Strategies
The model performance was evaluated using three strategies and was quantified using two statistical metrics including absolute root-mean-square errors (RMSE) and R 2 [51].
The ATC_GL combines both the spatially local information and LST-related descriptors based on the existing annual temperature cycle model. Three strategies were employed to evaluate the contributions of spatially local information, temporally global information, and LST-related descriptors, independently. Strategy 1 compares the performances of the ATC_GL and two classical ATC models (i.e., the ATCE and ATCT) when LSTs are randomly missing. Given that, satellite-derived LST images are invalid in a relatively Land 2023, 12, 309 9 of 25 large area. Strategy 2 compares the performances of the ATC_GL with classical temporally local interpolations, spatially local interpolations, and hybrid methods with varying LST gaps. Strategy 3 evaluated the model performances of the ATC_GL and several popular gap-filling methods under overcast conditions. Strategy 1: For the scenario with randomly missing LSTs under clear-sky conditions, the performance of the ATC_GL was compared with that of two typical ATC models (i.e., the ATCE and ATCO). For this strategy, 70% of the clear-sky LSTs were randomly selected as training data, and the rest were used for validation. This strategy considered data for Regions A-G.
Strategy 2: For the scenario with LST gaps of various sizes under clear-sky conditions, the performance of the ATC_GL was compared with that of several commonly used LST gap-filling methods, including the ATCH, Kriging interpolation, RSDAST, and HIT. The RSDAST, referring to a representative of spatiotemporally local methods, uses the distance relationships between invalid pixels and adjacent LSTs in spatiotemporally local windows to reconstruct missing LSTs [6]. The HIT is a hybrid method that integrates the observed values at the other overpasses in the same day, the spatiotemporally local interpolation, and the temporally global interpolation [18]. Ten LST images (each of a single day) of Region H, with more than 80% valid pixels, were randomly selected. We then produced 14 gaps with different sizes in the selected LST images, which are labelled as G01 (10 × 10) to G14 (500 × 500). The sizes of the produced gaps were 10 × 10, 20 × 20, 30 × 30, 40 × 40, 50 × 50, 60 × 60, 70 × 70, 80 × 80, 90 × 90, 100 × 100, 200 × 200, 300 × 300, 400 × 400, and 500 × 500 pixels. The performances of the methods above were assessed by the original valid LST values in the gaps.
Strategy 3: For the scenario of overcast conditions, Region A was selected to evaluate the performances of the ATC_GL and several popular LST gap-filling methods (including the ATCH, ATCE, Kriging, RSDAST, and HIT). For this strategy, the clear-sky values of LST MW in Region A were used as training data, while the LST MW values labelled with a cloud tag were used for validation. In addition, in situ measurements from seven SURFRAD sites in 2018 were used to evaluate the ATC_GL under overcast conditions.

Spatial Distribution of Model Performances
The overall performances of the ATC_GL, ATCE, and ATCO, as well as the accuracy improvement of the ATC_GL when compared with the ATCE (termed Drmse_ATCE) and ATCO (termed Drmse_ATCO), are shown in Figure 4. The mean RMSEs of the ATC_GL, ATCE, and ATCO models were 1.1 K, 3.5 K, and 4.2 K, respectively. As for different regions, the higher RMSEs of the ATC_GL occur in Regions B and E, and the relatively lower ones occur in the other regions. This may be because Regions B and E are located in the north and east of the Tibet Plateau, which is the junction of China's first and second ladders. The topography and climate conditions vary greatly in Regions B and E, causing drastic variations between adjacent pixels and days [12,51]. Consequently, the performance of the ATC_GL is lower in these high-heterogeneity regions than in other regions. Additionally, the RMSE distribution of the ATC_GL shows a smaller range (0.7-2.0 K) than that of the ATCE and ATCO (2.0-6.5 K; Figure 5) and the maximum frequency interval of the RMSE of each model is slightly smaller than the mean value of RMSE. The spatial distributions of the ATC_GL were relatively homogeneous within each study area when compared with those of the ATCE and ATCO. This result indicates that the ATC_GL can better predict spatiotemporally local variations during LST reconstruction than the other two ATC models. In Regions A, B, and C (located in northern, northeastern, and northwestern China, respectively), the accuracy of the ATC_GL was more than 3.0 K greater than that for the ATCE and ATCO ( Figure 4A4-C4,A5-C5); this result may be attributed to the high daily and intra-annual LST variation in these regions, where barren and grassland areas are the major land cover types [52][53][54][55].   The ATC_GL had the highest accuracy for all land cover types ( Figure 6). The highest and lowest RMSE for the ATC_GL was estimated in barren land (1.81 K) and farmland (0.92 K), respectively. The RMSEs of the ATCE and ATCO were relatively similar ( Figure  6); both the ATCE and ATCO showed maximum RMSEs in barren land (5.29 K and 5.53 K, respectively) and minimum RMSEs over water bodies (2.81 K and 3.04 K, respectively). The pronounced differences in RMSE between the land cover types for the ATCE and ATCO (more than 2.0 K) indicate that the performance of the ATC_GL remained relatively stable. The accuracy enhancement of the ATC_GL was the highest in barren areas (exceeding 3.5 K), with lower differences in accuracy observed between the models for forests, wetlands, and water bodies (more than 1.5 K; Figure 6). This result may reflect the reduced degree of LST fluctuation over dense vegetation and areas with high soil moisture because of higher levels of evapotranspiration and specific heat capacity [52][53][54][55]. The ATC_GL had the highest accuracy for all land cover types ( Figure 6). The highest and lowest RMSE for the ATC_GL was estimated in barren land (1.81 K) and farmland (0.92 K), respectively. The RMSEs of the ATCE and ATCO were relatively similar ( Figure 6); both the ATCE and ATCO showed maximum RMSEs in barren land (5.29 K and 5.53 K, respectively) and minimum RMSEs over water bodies (2.81 K and 3.04 K, respectively). The pronounced differences in RMSE between the land cover types for the ATCE and ATCO (more than 2.0 K) indicate that the performance of the ATC_GL remained relatively stable. The accuracy enhancement of the ATC_GL was the highest in barren areas (exceeding 3.5 K), with lower differences in accuracy observed between the models for forests, wetlands, and water bodies (more than 1.5 K; Figure 6). This result may reflect the reduced degree of LST fluctuation over dense vegetation and areas with high soil moisture because of higher levels of evapotranspiration and specific heat capacity [52][53][54][55]. The ATC_GL had the highest accuracy for all land cover types ( Table 4). The highest and lowest RMSE for the ATC_GL was estimated in barren land (1.8 K) and cropland (0.9 K), respectively. The RMSEs of the ATCE and ATCO were relatively similar (Table 4);  The ATC_GL had the highest accuracy for all land cover types ( Table 4). The highest and lowest RMSE for the ATC_GL was estimated in barren land (1.8 K) and cropland (0.9 K), respectively. The RMSEs of the ATCE and ATCO were relatively similar (Table 4); both the ATCE and ATCO showed maximum RMSEs in barren land (5.3 K and 5.9 K, respectively) and minimum RMSEs over water bodies (2.8 K and 3.3 K, respectively). The pronounced differences in RMSE between the land cover types for the ATCE and ATCO (more than 2.0 K) indicate that the performance of the ATC_GL (less than 0.9 K) remained relatively stable. The accuracy improvement of the ATC_GL was the highest in barren areas (exceeding 3.5 K), while there were relatively lower accuracies for forests, wetlands, and water bodies (more than 1.5 K; Table 4). The R 2 of the ATC_GL (exceed 0.94) was significantly higher than that of the ATCE (0.86-0.92) and the ATCO (0.87-0.94). A smaller RMSE usually corresponds to a larger R 2 for various land cover types. This result may reflect that the magnitude of the accuracy improvement depends on the LST fluctuation magnitude of the ATC_GL method. The improvement is relatively small where LSTs keep relatively stable, such as in forest, wetland and water body, probably because the ATCE and ATCO can achieve high accuracy in these areas [19]. On the contrary, the improvement is relatively higher where LSTs have high variations [56][57][58]. That is because the ATCE and ATCO have poor performance in such areas, while the ATC_GL can better reflect the spatial variations of LSTs in these areas by incorporating spatially local information. On a monthly scale, the performance of the ATC_GL was the most stable when compared with that of the ATCE and ATCO. Specifically, The RMSEs of the ATCE and ATCO were 2.6-4.4 K and 3.2-5.2 K, respectively, far exceeding the RMSE of the ATC_GL (0.9-1.7 K). Additionally, the RMSEs of the three models were relatively low from October to November and relatively high from May to July (Figure 7). The main reasons are as follows: the performances of the three models are generally related to the LST missing rate of the original MODIS LST data. There is no adequate neighborhood information when the LST missing rate is high; as a consequence, the spatiotemporal window will increase until the number of observed LSTs is satisfied. However, the use of LSTs with a large spatial or temporal distance can cause more uncertainties. Secondly, we speculate the performances of the three models are related to the LST values. The LSTs are high in the Northern Hemisphere in May to September, when the corresponding RMSEs of models are relatively large. When the LSTs are high, an annual temperature cycle may not be fully capable of capturing the seasonal variation of LSTs; as a consequence, the corresponding model RMSEs are relatively large. spatial or temporal distance can cause more uncertainties. Secondly, we speculate the performances of the three models are related to the LST values. The LSTs are high in the Northern Hemisphere in May to September, when the corresponding RMSEs of models are relatively large. When the LSTs are high, an annual temperature cycle may not be fully capable of capturing the seasonal variation of LSTs; as a consequence, the corresponding model RMSEs are relatively large. The ATC_GL can more accurately model detailed fluctuations within short intervals than the ATCE and ATCO, with the estimated values being largely consistent with the observed LSTs (Figure 8). The ATCO can effectively quantify intra-annual LST dynamics, but is less capable of handling daily LST fluctuations. The ATCE is superior to the ATCO in describing daily LST fluctuations (Figure 8b1-d1); however, the accuracy of the estimated fluctuations is affected by overestimation or underestimation. These results indicate that, although model performance can be enhanced by the inclusion of LST-related descriptors with temporally global interpolation, the incorporation of spatiotemporally global and local information along with appropriate statistical methods (i.e., RF) can more effectively improve model performance (refer to Section 5.1). The ATC_GL can more accurately model detailed fluctuations within short intervals than the ATCE and ATCO, with the estimated values being largely consistent with the observed LSTs (Figure 8). The ATCO can effectively quantify intra-annual LST dynamics, but is less capable of handling daily LST fluctuations. The ATCE is superior to the ATCO in describing daily LST fluctuations (Figure 8b1-d1); however, the accuracy of the estimated fluctuations is affected by overestimation or underestimation. These results indicate that, although model performance can be enhanced by the inclusion of LST-related descriptors with temporally global interpolation, the incorporation of spatiotemporally global and local information along with appropriate statistical methods (i.e., RF) can more effectively improve model performance (refer to Section 5.1).

Comparison of Model Performances in Strategy 2
The ATC_GL had the highest accuracy in most cases, and its performance was less affected by the size of the LST gap ( Figure 9). The RMSEs of the ATC_GL increased from 0.8 to 2.6 K with an increase in LST gap size. The performance of the Kriging interpolation decreased gradually with an increase in gap size, although it remained stable (1.4-2.0 K) at relatively small gap sizes (smaller than 50 × 50 pixels). The RMSEs of the HIT were 1.0-1.5 K higher than those of the ATC_GL, although the variation across time was similar between the models. The accuracy of the RSDAST was slightly higher than that of the ATC_GL at very small gap sizes (smaller than 20 × 20 pixels). However, its performance decreased dramatically with an increase in gap size, and the RSDAST was the poorest performing model at gap sizes larger than 80 × 80 pixels.  (a-d)). Each randomly selected pixel corresponds to a specific subfigure. Subfigures (a1-d1) are the enlarged insets of subfigures (a-d).

Comparison of Model Performances in Strategy 2
The ATC_GL had the highest accuracy in most cases, and its performance was less affected by the size of the LST gap ( Figure 9). The RMSEs of the ATC_GL increased from 0.8 to 2.6 K with an increase in LST gap size. The performance of the Kriging interpolation decreased gradually with an increase in gap size, although it remained stable (1.4-2.0 K) at relatively small gap sizes (smaller than 50 × 50 pixels). The RMSEs of the HIT were 1.0-1.5 K higher than those of the ATC_GL, although the variation across time was similar between the models. The accuracy of the RSDAST was slightly higher than that of the ATC_GL at very small gap sizes (smaller than 20 × 20 pixels). However, its performance Compared to the other models, the ATC_GL had the highest accuracy (RMSE = 1.7 K) and preserved the best spatial detail at the largest LST gap size (500 × 500 pixels; Figure 10). The RSDAST (RMSE = 3.0 K) and Kriging interpolation (RMSE = 3.4 K) showed the worst performance in estimating spatial variation. These findings indicate that spatially local interpolations may be less effective or even unfeasible at large gap sizes. The HIT is relatively capable of estimating the LST contrast in this region; however, its overall accuracy was low (RMSE = 2.9 K). The ATCH (RMSE = 2.0 K) had the second highest accuracy among the models and showed the greatest differences in accuracy with the ATC_GL in the regions as shown by the red boxes in Figure 10g.
decreased dramatically with an increase in gap size, and the RSDAST was the poorest performing model at gap sizes larger than 80 × 80 pixels. Compared to the other models, the ATC_GL had the highest accuracy (RMSE = 1.7 K) and preserved the best spatial detail at the largest LST gap size (500 × 500 pixels; Figure  10). The RSDAST (RMSE = 3.0 K) and Kriging interpolation (RMSE = 3.4 K) showed the worst performance in estimating spatial variation. These findings indicate that spatially local interpolations may be less effective or even unfeasible at large gap sizes. The HIT is relatively capable of estimating the LST contrast in this region; however, its overall accuracy was low (RMSE = 2.9 K). The ATCH (RMSE = 2.0 K) had the second highest accuracy among the models and showed the greatest differences in accuracy with the ATC_GL in the regions as shown by the red boxes in Figure 10g.

Comparison of Model Performances in Strategy 3
For the scenario of overcast conditions, the performance of the ATC_GL was again compared with those of the LST gap-filling methods mentioned in Strategy 2 (the ATCE, ATCH, RSDAST, Kriging, and HIT) (refer to Appendix C for the details of study areas and data). The ATC_GL had the highest accuracy in overcast conditions (RMSE = 3.7 K, R 2 =

Comparison of Model Performances in Strategy 3
For the scenario of overcast conditions, the performance of the ATC_GL was again compared with those of the LST gap-filling methods mentioned in Strategy 2 (the ATCE, ATCH, RSDAST, Kriging, and HIT) (refer to Appendix C for the details of study areas and data). The ATC_GL had the highest accuracy in overcast conditions (RMSE = 3.7 K, R 2 = 0.93; Figure 11). The RMSEs of the ATCE, ATCH, Kriging interpolation, and HIT were 5.8 K, 5.1 K, 4.7 K, and 5.2 K, respectively. Notably, the LSTs in overcast conditions may be overestimated by the ATC_GL especially in the range of 275.0 K to 290.0 K, although the level of overestimation is relatively lower when compared with other gap-filling models ( Figure 10). The overestimation may be associated with the inadequate capability of SAT to fully grasp the LST dynamics under overcast conditions, especially for regions with large LST gaps [19,59]. We further need to note that the ATC_GL is indeed a framework for LST reconstruction, which can adjust parameter combinations in different scenarios. The performance of the framework is expected to improve if more auxiliary parameters related to overcast conditions (e.g., longwave and shortwave radiation and soil moisture) are incorporated [41]. In addition, the performance of the ATC_GL was evaluated using SURFRAD LSTs under both clear-sky and overcast conditions. Prior to the model evaluation, the systematic differences between MODIS and SURFRAD LSTs were corrected mainly by referring to previous studies [35,60,61]. The ATC_GL maintains high accuracy under both clear-sky (RMSE ranges from 0.5 to 1.0 K) and overcast conditions (RMSE ranges from 2.1 to 4.1 K) Figure 11. Accuracy evaluations of the reconstructed LSTs using different methods under overcast conditions. Subfigure (g) shows the bias of different methods in different LST ranges.
In addition, the performance of the ATC_GL was evaluated using SURFRAD LSTs under both clear-sky and overcast conditions. Prior to the model evaluation, the systematic differences between MODIS and SURFRAD LSTs were corrected mainly by referring to previous studies [35,60,61]. The ATC_GL maintains high accuracy under both clear-sky (RMSE ranges from 0.5 to 1.0 K) and overcast conditions (RMSE ranges from 2.1 to 4.1 K) ( Table 5). We believe that the performance of the ATC_GL may be further improved once more auxiliary data under overcast conditions are employed.

Advantages of the ATC_GL Reference to Previous Methods
Our assessments have shown that the ATC_GL has the highest accuracy when compared with several commonly used gap-filling methods (i.e., the Kriging interpolation, RSDAST, ATCE, ATCH, and HIT). The RSDAST and Kriging methods mainly use spatially or spatiotemporally local information for LST reconstruction and exhibit relatively high accuracy when LST gaps are small. However, as LST gap size increases, the performance of these two methods decreases significantly, likely due to the lack of neighborhood information and LST-related descriptors (e.g., NDVI and SAT), both of which are important factors in LST reconstruction. By comparison, the ATCE and ATCH do incorporate LST-related descriptors in a temporally global fashion, and their accuracies are generally insensitive to the size of the LST gap [19]. However, the ATCE and ATCH ignore spatiotemporally local information among adjacent pixels and only use linear regressions to characterize the relationships between LST and LST-related descriptors; thus, their accuracy is usually lower than those that incorporate spatially local information especially when the LST gaps are relatively small [42]. For LST reconstruction, the HIT incorporates a dynamic spatiotemporal window to avoid error accumulation during repeated interpolations with the already reconstructed LSTs [18]. However, similar to the RSDAST and Kriging methods, the HIT does not use LST-related descriptors and therefore has a relatively low performance under overcast conditions. The ATC_GL incorporates both temporally global and spatiotemporally local interpolations. From the temporally global perspective (i.e., monthly or seasonal cycles), the ATC_GL incorporates the seasonal cycle of solar radiation to reconstruct intra-annual LST dynamics; therefore, it is insensitive to LST gap size, especially compared with the methods based on spatiotemporally local interpolations. From the spatiotemporally local scale, the ATC_GL incorporates the RF algorithm in obtaining the complex relationships between daily LST fluctuations and a series of relevant factors (e.g., the daily SAT fluctuations, NDVI, and DEM) within dynamic spatiotemporal windows. The ATC_GL is therefore suitable for LST reconstruction with complex surface properties. Briefly, the ATC_GL presents the following two key advantages when compared with previous ATC models: (1) the better use of spatiotemporally local information within dynamic spatiotemporal windows, and (2) a more realistic estimation of the relationships between LSTs and LST-related descriptors with the RF algorithm. The ATC_GL is a simple yet robust model that is able to reconstruct accurately spatiotemporally continuous LSTs. In addition, the ATC_GL is also independent of auxiliary data (i.e., only requiring globally acquirable LST and reanalysis data) and thus applicable globally.

Contributions of LST-Related Descriptors in the Estimation of Daily LST Fluctuations
We quantified the contributions of the LST-related descriptors used in ATCT (i.e., daily SAT fluctuation, NDVI, and DEM) with an attribution approach [62]. The daily SAT fluctuation had the highest mean annual contribution (40.6%), followed by NDVI (34.9%) and DEM (24.5%; Figure 12). These results also highlight the deterministic role of daily SAT fluctuation and the indispensable contributions of NDVI and DEM in LST reconstruction. The relevance of these factors to LST reconstruction comes from: (1) the close correlation between LST and SAT under different conditions [63,64], and (2) the regulatory effects of NDVI and DEM on the daily variation of LST [59,63]. Furthermore, notable monthly variations were observed in the contributions of these three variables: (1) The contribution of daily SAT fluctuation was relatively high from May to September (exceeding 40.0%). This was likely owing to the increased amount of cloud coverage during this period that consequently weakened the regulation of LST variation by NDVI and DEM. (2) The contributions of NDVI were relatively higher in January, April, October, and December than those in the other months; this finding could be explained by sudden changes in plant cover (e.g., crop harvesting and sowing). (3) The contribution of the DEM was relatively stable within the annual cycle, with the largest inter-month difference of lower than 15.0%, probably because DEM remained unchanged throughout the year.

Limitations and Prospects
Regardless of its comparatively high performance in the modeling of intra-annual LST dynamics, the ATC_GL has some limitations. First, the ATC_GL ignores the overpassing time difference of MODIS LSTs (up to 2 h) within a daily cycle [65] and viewing angle differences among temporally adjacent days [66][67][68]. Although previous studies have shown that the accuracy enhancement after the normalization of overpassing time and angle is limited [19], the ATC_GL is still expected to benefit if these factors are considered.
Second, the ATC_GL is prone to overestimate LSTs in overcast conditions, probably because the used LST-related descriptors have less prediction ability in such conditions. Microwave data can be used to directly obtain the surface thermal status under overcast conditions and, therefore, could be useful for LST reconstruction. In situ measurements

Limitations and Prospects
Regardless of its comparatively high performance in the modeling of intra-annual LST dynamics, the ATC_GL has some limitations. First, the ATC_GL ignores the overpassing time difference of MODIS LSTs (up to 2 h) within a daily cycle [65] and viewing angle differences among temporally adjacent days [66][67][68]. Although previous studies have shown that the accuracy enhancement after the normalization of overpassing time and angle is limited [19], the ATC_GL is still expected to benefit if these factors are considered.
Second, the ATC_GL is prone to overestimate LSTs in overcast conditions, probably because the used LST-related descriptors have less prediction ability in such conditions. Microwave data can be used to directly obtain the surface thermal status under overcast conditions and, therefore, could be useful for LST reconstruction. In situ measurements under overcast conditions, such as the LST or other remote sensing data that are directly related to LST (e.g., radiation), could also be helpful to minimize the overestimation [10,[13][14][15]41,[69][70][71][72][73][74][75]. The performance of the ATC_GL is expected to improve once these auxiliary data are integrated. Fortunately, the ATC_GL is not a specific algorithm but a framework that can incorporate both spatiotemporally global and local information, for which auxiliary data can be integrated directly as LST-related descriptors during LST reconstruction. When applying the ATC_GL to LST reconstruction, practitioners should consider both the availability of LST-related descriptors and the model efficiency. The ATC_GL maintains high accuracy under both clear-sky (RMSE ranges from 0.5 to 1.0 K) and overcast conditions (RMSE ranges from 2.1 to 4.1 K). We believe that the performance of the ATC_GL may be further improved once more auxiliary data under overcast conditions are employed. In addition, the ATC_GL is expected to be applicable for reconstructing LSTs with long time gaps such as the Landsat thermal data [28,70] (refer to Appendix C for the details).

Conclusions
The existing ATC models are temporally global interpolations of nature, without considering the spatiotemporally local information among adjacent pixels as well as the complex relationship between LST and LST-related descriptors. These limitations translate to relatively low-performance estimates under varying conditions. To overcome these limitations, we proposed an improved ATC model (i.e., the ATC_GL) that integrates both spatiotemporally global and local information and incorporates the RF algorithm and a series of auxiliary data (e.g., daily SAT variation and NDVI) to model the complex relationships between LSTs and LST-related descriptors.
Our assessments showed that, in the scenario with randomly missing LSTs under clear-sky conditions, the RMSE of the ATC_GL was 1.1 K (2.3 K and 3.1 K lower than that of the ATCE and ATCO, respectively). The performance of the ATC_GL was more stable between land cover types when compared with the ATCE and ATCO. For the scenario with LST gaps of various sizes under clear-sky conditions, our assessments revealed that the accuracy of the ATC_GL was generally the highest when compared with several LST gap-filling methods. The ATC_GL is a framework rather than a specific algorithm, as it can directly integrate relevant LST-related descriptors in LST reconstruction. We believe that the ATC_ GL would improve the accuracy of intra-annual LST dynamics estimation, LST reconstruction, and other related applications.

Appendix A
The size of the initial spatiotemporal window can affect model accuracy and efficiency [6,71,72]. In order to identify the optimal spatiotemporal window, we examine the accuracy variations of the ATC_GL depending on the increase in spatial window size (increasing from 5 × 5 pixels to 19 × 19 pixels, Figure A1) of Region H and the temporal window size (increasing from 1 to 9 days, Figure A2). The assessments reveal that the accuracy reaches the highest when the initial spatial window sizes are 9 × 9 pixels and 11× 11 pixels (i.e., w3 and w4 as given in Figure A1); and the associated mean RMSE for these two window sizes is around 1.3 K. We therefore chose the size of the initial spatiotemporal window as 9 × 9 pixels for the ATC_GL. The model achieves the best performance when the initial temporal window is 1 day with a mean RMSE of 1.2 K.  Figure A1. The accuracy of the ATC_GL along with the increase in the size of the initial spatial window. w1-w8 denote the size of 5 × 5, 7 × 7, 9 × 9, 11 × 11, 13 × 13, 15 × 15, 17 × 17, and 19 × 19 pixels of initial spatial window, respectively. Figure A2. The accuracy variations of the ATC_GL depending on the increase in the size of the initial temporal window. T1-T8 denote the size of 1-9 days for the initial temporal window, respectively.
The ATC_GL performance shows different sensitivity in response to temporal and spatial window sizes (Figures A1 and A2). Therefore, the different combinations of increased spatial and temporal window sizes may affect the ATC_GL performance. To in- Figure A1. The accuracy of the ATC_GL along with the increase in the size of the initial spatial window. w1-w8 denote the size of 5 × 5, 7 × 7, 9 × 9, 11 × 11, 13 × 13, 15 × 15, 17 × 17, and 19 × 19 pixels of initial spatial window, respectively.  Figure A1. The accuracy of the ATC_GL along with the increase in the size of the initial spatial window. w1-w8 denote the size of 5 × 5, 7 × 7, 9 × 9, 11 × 11, 13 × 13, 15 × 15, 17 × 17, and 19 × 19 pixels of initial spatial window, respectively. Figure A2. The accuracy variations of the ATC_GL depending on the increase in the size of the initial temporal window. T1-T8 denote the size of 1-9 days for the initial temporal window, respectively.
The ATC_GL performance shows different sensitivity in response to temporal and spatial window sizes (Figures A1 and A2). Therefore, the different combinations of increased spatial and temporal window sizes may affect the ATC_GL performance. To investigate the impact of different combinations of increased spatial and temporal window sizes on the ATC_GL performance, a new experiment was conducted. The results show Figure A2. The accuracy variations of the ATC_GL depending on the increase in the size of the initial temporal window. T1-T8 denote the size of 1-9 days for the initial temporal window, respectively. The ATC_GL performance shows different sensitivity in response to temporal and spatial window sizes (Figures A1 and A2). Therefore, the different combinations of increased spatial and temporal window sizes may affect the ATC_GL performance. To investigate the impact of different combinations of increased spatial and temporal window sizes on the ATC_GL performance, a new experiment was conducted. The results show that the accuracy of the ATC_GL rarely changes with different combinations of increased spatial and temporal window sizes (Table A1). When using a different combination of increased spatial and temporal window sizes, the ATC_GL performance varies by less than 14% when referenced to the original spatiotemporal window sizes (i.e., 2 km and 2 days). This suggests that the responsiveness of the ATC_GL performance to spatial and temporal window sizes would not largely bias the major conclusions. spatial and temporal window sizes (Table A1). When using a different combination of increased spatial and temporal window sizes, the ATC_GL performance varies by less than 14% when referenced to the original spatiotemporal window sizes (i.e., 2 km and 2 days). This suggests that the responsiveness of the ATC_GL performance to spatial and temporal window sizes would not largely bias the major conclusions.

Appendix C
Land surface temperatures may be missing for a long time interval, e.g., the LST missing gaps of the Landsat thermal data can be as large as half a month or more due to the long revisit period of Landsat data and cloud contamination. To test the performance of the ATC_GL when the LST missing gap is large, Region A was chosen as the study area for evaluating model performances when the LST missing gap is as large as one month. Our results show that the ATC_GL maintains a high accuracy when the LSTs are missing for a long time interval. For example, the tests using Region A as the study area by setting the LST missing time interval of one month indicate that the mean RMSE of the ATC_GL can reach 3.3 K (Figure A4), an accuracy that is generally acceptable especially when compared with other popularly used methods for such purposes. Therefore, the ATC_GL is expected to be applicable for reconstructing LSTs with long time gaps such as the Landsat thermal data. The ATC_GL should also be helpful for the generation of spatiotemporally seamless LST products on a daily basis at the spatial resolution of Landsat thermal data (i.e., 100 m).

Appendix C
Land surface temperatures may be missing for a long time interval, e.g., the LST missing gaps of the Landsat thermal data can be as large as half a month or more due to the long revisit period of Landsat data and cloud contamination. To test the performance of the ATC_GL when the LST missing gap is large, Region A was chosen as the study area for evaluating model performances when the LST missing gap is as large as one month. Our results show that the ATC_GL maintains a high accuracy when the LSTs are missing for a long time interval. For example, the tests using Region A as the study area by setting the LST missing time interval of one month indicate that the mean RMSE of the ATC_GL can reach 3.3 K (Figure A4), an accuracy that is generally acceptable especially when compared with other popularly used methods for such purposes. Therefore, the ATC_GL is expected to be applicable for reconstructing LSTs with long time gaps such as the Landsat thermal data. The ATC_GL should also be helpful for the generation of spatiotemporally seamless LST products on a daily basis at the spatial resolution of Landsat thermal data (i.e., 100 m).