Quantifying the Inﬂuences of Driving Factors on Land Surface Temperature during 2003–2018 in China Using Convergent Cross Mapping Method

: The relationship between land surface temperature (LST) and environmental factors is complex and nonlinear. To determine these relationships for China, this study analyzed the driving effects of air temperature, vegetation index, soil moisture, net surface radiation, precipitation, aerosols, evapotranspiration, and water vapor on LST based on remote-sensing and reanalysis data from 2003–2018, using a convergent cross-mapping method. During the study period, air temperature and net surface radiation were the dominant drivers of LST with a cross-mapping skill above 0.9. Vegetation index and evapotranspiration were the secondary drivers of LST with a cross-mapping skill that was higher than 0.5. Except for air temperature and net surface radiation, the direction and strength of the effects of the driving factors on LST were related to the climate type. The effects of air temperature and net radiation on LST diminished from north to south, indicating that LST was more sensitive to air temperature and net radiation in energy-limited regions. However, the effects of vegetation index and evapotranspiration on LST varied signiﬁcantly across climate zones; that is, positive effects were mostly in non-monsoonal zones and negative effects were primarily in monsoonal zones. Our results quantiﬁed the driving role of environmental factors on LST and provided a comprehensive understanding of LST dynamics.


Introduction
In a complex climate system, anomalies in key variables often lead to chain effects, disrupt related variables, and even cause fluctuations in the entire system. Land surface temperature (LST), a central element in the climate system, plays a fundamental role in assessing land-atmosphere interactions and climate change. It is defined as a high priority parameter by the International Geosphere-Biosphere Programme and as an important climate variable by the Global Climate Observing System (GCOS) of the World Meteorological Organization [1]. LST anomalies have a significant impact on energy cycles, ecological balance, and even human life. Unfortunately, the drivers that cause changes in LST are not accurately determined. Therefore, detecting the causal relationship between drivers and LST is crucial for analyzing both temperature anomalies and climate change.
Over the past few decades, significant results have been obtained for LST and its drivers, such as the relationship between surface cover type, NDVI (Normalized Difference Vegetation Index), elevation, soil moisture, evapotranspiration, and surface radiation with LST [2][3][4][5][6]. Generally, studies have focused on the urban heat island effect, heat wave phenomenon, or the effect of vegetation changes on LST [7][8][9][10]. There are limited studies analyzing the relationship between LST and environmental factors from causality perspective, especially the analysis of multiple environmental factors and LST. Several non-causal methods have been frequently used to analyze the relationship between environmental factors and LST in previous studies, including the correlation coefficient method [11], linear regression [12], stepwise regression [13,14], and geographically weighted regression [15,16]. In addition, most studies have statistically analyzed the correlations between LST and environmental factors on a time series, assuming that the environmental factors with significant correlations are the drivers of LST. However, they have ignored the basic principle that correlations are not causal. Therefore, when studying the interaction between LST and environmental factors, it is important to determine whether there is a causal relationship between the two.
In recent years, there has been a proliferation of approaches to causality analysis, which can be classified into qualitative process models and data-driven models. Qualitative process models require sufficient theory and experience to analyze the causality of things; however, obtaining accurate and sufficient process knowledge is difficult and time-consuming. By contrast, data-driven approaches can provide multiple ways to find causal relationships based on large amounts of historical data. Data-driven methods have been widely applied to identify the causal relationships between variables [17][18][19][20][21][22][23][24][25]. The Granger causality (GC) test proposed by Granger [26] is considered to be one of the earliest data-driven causality detection methods. The GC test assumes that a causal relationship exists between two variables if the prediction of one variable can be significantly improved by combining the information from another variable. However, GC tests require variables to be independent of each other; therefore, they are not suitable for causality detection of nonlinear systems with coupled relationships. The geographic detector method proposed by Wang et al. [27] determines the role of influencing factors through the spatial divergence between the independent and dependent variables, which not only detects the influencing factors, but also provides the magnitude of the ability of the influencing factors to explain the spatial divergence of the dependent variables. However, it has certain requirements with respect to the data types of variables, especially the causality of continuous time series variables detection. Thus, the results can be anomalous.
The empirical dynamic modeling (EDM) causality detection method was proposed by Sugihara et al. [17]. This convergent cross mapping (CCM) method is a practical databased causality detection method based on Takens' theorem to identify causal relationships in complex nonlinear systems with coupling. According to historical studies, the CCM method can effectively identify causal relationships in non-separable nonlinear systems, while being an equation-free statistical method that is not controlled by models [28]. The CCM method can not only identify causal relationships between variables while effectively eliminating the effects of other factors, but can also eliminate pseudo-causal relationships caused by synchronization effects [29]. In addition, the high sensitivity of CCM methods to weakly coupled relationships in nonlinear systems can lead to their increasing interest and applications in various fields, such as fisheries, medicine, chemistry, climate change, environmental pollution, and ecology [30][31][32][33][34]. However, to the best of our knowledge, few studies have focused on the causal relationship between LST and environmental elements based on the CCM method.
The main objectives of this study were to (1) identify the causal relationships between LST and environmental factors based on the CCM method, (2) screen and eliminate pseudocausal relationships caused by synchronous effects between LST and environmental factors, and (3) explore the direction and intensity of the driving factors on LST using scenario exploration. The remainder of this paper is organized as follows. Section 2 presents the data sources for the environmental factors. Section 3 introduces the methodology. Section 4 provides the results before the discussion in Section 5, and Section 6 concludes the study. The MYD11C3 product was obtained by synthesizing and averaging the values of the corresponding months in the MYD11C1 product. This product provided LST estimates for day (1:30 p.m.) and night (1:30 a.m.) overpass times [35][36][37][38]. This study limited the analysis to high-quality clear-sky pixels with an LST error ≤ 1 K using the MYD11C3 quality control (QC) layer. The MODIS Aqua overpass time roughly corresponded to the daily maximum LST, which occurred in the afternoon, and the daily minimum LST, which was after midnight. Therefore, the average of the two LSTs was used in this study.
As an indicator closely related to green biomass and leaf area indices at regional and global scales, NDVI has been widely recognized for studying the relationship between LST and vegetation [39]. The MYD13C2 product is analogous to MYD13C1 (Vegetation Indices 16-Day L3 Global 0.05 • CMG) but is based on MYD13A3 (Vegetation Indices Monthly L3 Global 1 km SIN Grid) for a monthly temporal resolution. All other specs are maintained and the production features are retained. In this study, the MYD04_L2 and MYD05_L2 products were used to obtain AOD (Aerosol Optical Depth) and WV (Water vapour) factors. The selected bands of the two products were AOD_550_Dark_Target_Deep_Blue_Combined and water vapor infrared.
The MYD16A2 product was an 8-day composite dataset produced at a spatial resolution of 500 m. The algorithm used for the MYD16 product was based on the logic of the Penman-Monteith equation, which included inputs of daily meteorological reanalysis data along with Moderate Resolution Imaging Spectroradiometer (MODIS) remotely sensed data products, such as vegetation property dynamics, albedo, and land cover [40].

ESA-CCI Soil Moisture Data
The soil moisture dataset was provided by the European Space Agency (ESA) (http://esa-soilmoisture-cci.org (accessed on 1 November 2021)). ESA released the global soil moisture product, which included active microwave products, passive microwave products, and combined active-passive microwave products [41]. In this study, the combined active-passive microwave product was analyzed from 2003 to 2018. The dataset had a spatial resolution of 0.25 • × 0.25 • [42]. Moreover, its applicability in China was demonstrated by An et al. (2016) [43].

ERA5-Land and CRU Data
ERA5-Land is a reanalysis dataset that provides a consistent view of the evolution of land variables over several decades at an enhanced resolution compared with ERA5. ERA5-Land was produced by replaying the land component of the European Center for Medium-Range Weather Forecast (ECMWF) ERA5 climate reanalysis [44]. Three products were used in this study, including air temperature at 2 m above the ground, surface net solar radiation, and surface net thermal radiation, with a spatial resolution of 0.1 • × 0.1 • . The surface net solar radiation and surface net thermal radiation data were used to calculate the surface net radiation.
Monthly CRU dataset was found to have a better temporally continuous availability than traditional weather station observations [45].
All data from 2003 to 2018 were spatially and temporally scaled at a uniform resolution of 0.05 • and monthly time steps, respectively. Spatial resampling to a 0.05 • resolution was performed by applying bilinear interpolation. Daily and 8-day composite data were averaged to obtain the monthly values. Data with zero values are treated as invalid pixels. A summary of these data is presented in Table 1.

Method
To analyze the causality between LST and its related environmental factors, five main procedures for causality detection were summarized as follows: (1) The CCM method was used to identify the causality between LST and environmental factors. (2) The seasonal surrogate test (SST) was used to remove pseudo-causality between the LST and environmental factors caused by the seasonal synchronization effect. (3) Multivariate EDM helped examine whether environmental factors had a stronger causal effect on LST. (4) The ECCM (Extend Convergence Cross Mapping) method was used to further identify pseudo-causality due to generalized synchronization and Moran effects. (5) Multivariate scenario exploration was employed to determine the direction and strength of the real driver on LST. Details of each procedure are described in the following sections. A flowchart of causality detection between the LST and environmental factors is shown in Figure 1.

CCM
The CCM provides a way to identify the causal relationship between two variables in a complex nonlinear system. The fundamental principle of the CCM method is as follows.
If variable x in a system acts on variable y, the historical information of x is recorded by y. The time series of y can be used to recover the historical information of x; that is, the historical values of y can be used to map the time series of x. Details can be found in the three animations provided by Sugihara et al. [17] (http://science.sciencemag.org/ content/suppl/2012/09/19/science.1227079.DC1 (accessed on 1 November 2021)). For the CCM analysis, the simplex projection was used to predict x with the historical information of y [46]. Cross-mapping skill (ρ ccm ) was defined as the Pearson correlation coefficient ρ between the predicted and actual values of x. If the value of ρ ccm increased with the length of the time series and convergence was achieved, the causal effect of x on y could be determined. The Mann-Kendall monotonicity detection method (M-K test) was used to determine the convergence of ρ ccm [47,48]. The validity of the CCM detection results is determined by the convergence of ρccm, whereby the CCM method requires the time series length of the data to be greater than 30 [17]. To further improve the accuracy of the results in this study, we only detect the pixels with time series length greater than 90.
To explore the driving factors of LST, eight environmental factors including air temperature (TA), normalized difference vegetation index (NDVI), soil moisture (SM), net surface radiation (RN), precipitation (PRE), aerosol optical depth (AOD), evapotranspiration (ET), and water vapor (WV) were selected as potential driving variables. In this study, the CCM method was used to analyze whether the aforementioned eight environmental factors affected the response variable LST.

CCM
The CCM provides a way to identify the causal relationship between two variables in a complex nonlinear system. The fundamental principle of the CCM method is as follows.
If variable x in a system acts on variable y, the historical information of x is recorded by y. The time series of y can be used to recover the historical information of x; that is, the historical values of y can be used to map the time series of x. Details can be found in the three animations provided by Sugihara et al. [17] (http://science.sciencemag.org/content/suppl/2012/09/19/science.1227079.DC1 (accessed on 1 November 2021)). For the CCM analysis, the simplex projection was used to predict x with the historical information of y [46]. Cross-mapping skill (ρccm) was defined as the Pearson correlation coefficient ρ between the predicted and actual values of x. If the value of ρccm increased with the length of the time series and convergence was achieved, the causal effect of x on y could be determined. The Mann-Kendall monotonicity detection method (M-K test) was used to determine the convergence of ρccm [47,48]. The validity of the CCM detection results is determined by the convergence of ρccm, whereby the CCM method requires the time series length of the data to be greater than 30 [17]. To further improve the accuracy of the results in this study, we only detect the pixels with time series length greater than 90.

SST
A similar seasonal cycle between LST and an environmental factor may result in a high value of ρ ccm , even if there is no causal relationship. This phenomenon is described as the seasonal synchronization effect [17,29]. In this case, the CCM method could not fully distinguish real causality or spurious correlation [49]. A null test with a surrogate time series was considered as an effective method to address this problem [20].
The time series of an environmental factor with seasonal cycle changes {Z} was decomposed into seasonal cycle series {Zs} and short-term fluctuation series {∆d}. {∆d} was randomly shuffled to obtain a new fluctuation series {∆d*}. Then, {∆d*} was added to the seasonal cycle series {Zs} to obtain the surrogate time series of the environmental factor, that is, {Z*} = {Zs} + {∆d*}. It has been established that {Z*} has the same seasonal cycle series as {Z} but a different fluctuation series. If {Z} is indeed the driver of LST, LST will be sensitive not only to the seasonal cycle series {Zs}, but also to the fluctuation series {∆d}. Therefore, LST should be better at predicting the actual series {Z} than the surrogate series {Z*}. The shuffling procedure was repeated 500 times to produce a set of surrogate series of environmental factors. The value of ρ ccm between the LST and each surrogate series was calculated. The environmental factor was retained if the value of ρ ccm between the LST Remote Sens. 2022, 14, 3280 6 of 26 and the actual series of environmental factors was significantly higher than that between the LST and the 500 surrogate series. Otherwise, the environmental factor was considered to have no causal relationship with LST and was rejected. Finally, a non-parametric onesample Wilcoxon test was chosen to estimate whether the actual series of environmental factors had a stronger effect on LST than the surrogate series at a significant level (p < 0.05).

Multivariate EDM
For causality detection among random variables, the CCM method may not be able to accurately identify the driving effect of an environmental factor on LST, due to data limitations, such as time series length and breakpoints in the time series [29]. In this study, multivariate EDM was used to determine whether an environmental factor had a significant driving effect on LST by comparing the improvement in EDM forecasting skills [50]. In other words, if an environmental factor had a driving effect on LST, it would significantly improve the forecasting skill of the model by adding environmental factors into the multivariate EDM model.
In this study, a multivariate EDM model composed of the time series of an environmental factor and the lagged time series of LST was used to forecast the LST. Forecasting skill (ρ) is the correlation coefficient between the time series of the forecasted and the original LST. Forecast improvement is the difference between the forecasting skill of the model with the environmental factor (ρ with_factor ) and that without the environmental factor (ρ without_factor ), that is, forecast improvement ∆ρ = ρ with_facto r − ρ without_facto r. A non-parametric one-sample Wilcoxon test was used to test the significance (p < 0.05) of the forecasting improvement of the multivariate EDM model.

ECCM
In a complex system, the accuracy of causality detection is affected by the generalized synchronization and the Moran effects [29], which cannot be removed by the SST and multivariate EDM methods. Ye et al. (2015) addressed the above problems by considering different time lags. In an actual causal relationship, the driving variable x can only act on the present or future values of the response variable y and not on the past values of y. Therefore, in the cross-mapping process, the response variable y can only predict the present or past values of the driving variable x, but not the future values of x. Accordingly, when using the CCM method to detect causality between variables x and y, if there is a time lag in the driving effect, the time lag must be non-positive. Based on this concept, Ye et al. [29] proposed an extension of the CCM method (ECCM). First, the set of lagged time series of variable x was established; ρ ccm between each lagged time series and variable y was calculated separately, and the time lag corresponding to the optimal ρ ccm was selected. If time lag ≤ 0, there was a driving effect of variable x on y, and the driving effect has a time lag of lag length. If time lag > 0, there is no driving effect of x on y. In this study, we performed ECCM detection for all environmental factors and LST.

Multivariate Scenario Exploration
After removing the pseudo-causal relationship between the environmental factors and LST, a scenario analysis was conducted to further explore the direction and strength of the driving effect of the true driver on LST. Multivariate scenario analysis based on EDM is a method for empirically assessing the effect of small changes in the driver on LST, which can effectively assess the response of LST to fluctuations in the driver [20]. The specific steps of multivariate scenario analysis are detailed below. Taking air temperature (TA) as an example, we predicted the effect of small changes in TA on LST for a corresponding lag length. The multivariate EDM model in this study included the lagged coordinates of the LST and the time series of TA: (1) This equation was obtained with S-maps, where E is the optimal number of coordinates for the model to make predictions. For each historical time point t, this study predicted the LST with a small increase (+∆TA/2) and a small decrease (−∆TA/2) in the historically measured driver TA(t). To compare the strength of the effect of TA on LST in different regions, ∆TA was used, which corresponded to~5% of the standard deviation of TA across China.
The value of ∆LST(t)/∆TA indicates the magnitude of the effect of TA on LST and the sign indicates the direction of the effect. The spatial distribution of LST showed an obvious latitudinal characteristic; that is, it decreased with increasing latitude. In addition, it showed different patterns influenced by altitude. For example, the Qinghai-Tibet Plateau in the southwest, which is the highest-altitude region in China, showed low temperatures at lower latitudes. The spatial distribution of TA was basically the same as that of LST, but it was approximately 5-9 K lower in magnitude than LST. NDVI and SM showed a distribution pattern that increased gradually from northwest to southeast, whereas RN and PRE showed a gradual increase from south to north. AOD was significantly higher in the northwest desert region and in the eastern region with higher population density; however, it was lower in all other regions. The magnitude of ET was strongly related to vegetation type; therefore, it was similar to NDVI in terms of spatial distribution. However, there was a large amount of missing data (blank areas in Figure 2h) in the northwest because the region was covered with deserts, especially Gobi. WV decreased from the coast to the inland and was at its lowest in the Qinghai-Tibet Plateau region, due to its leading altitude and topography.

Causality Detection Process between LST and Environmental Factors
To clearly express the process of causality detection between LST and environmental factors, a pixel with grassland land cover type from the PMC region was selected as an example for a detailed demonstration. To ensure more stable results, each pair of crossmapping processes was repeated 100 times [19], and the average of the 100 results was used as the final result. Figure 3 shows the detection results between LST and each environmental factor, with the horizontal axis being the length of the time series and the vertical axis being the ρ ccm between LST and environmental factors.
As shown in Figure 3, the ρ ccm curves between LST and each environmental factor exhibited different trends with an increase in the length of the time series. Among them, ρ ccm between RN and LST, TA and LST, and NDVI and LST increased rapidly with increasing time series length, and converged to 0.99, 0.97, and 0.91, respectively. This indicated strong potential causal relationships between RN and LST, TA and LST, and NDVI and LST. Although the ρ ccm between ET and LST, PRE and LST, and WV and LST showed monotonically increasing and converging characteristics, the corresponding ρ ccm values at convergence were relatively small, that is, 0.48, 0.42, and 0.38, respectively. This revealed that the strength of the potential causal relationships between them and LST was relatively low. The ρ ccm values between SM and LST, and AOD and LST decreased with increasing time series, indicating that there was no causal relationship between them and LST. The monotonic convergence of ρ ccm between environmental factors and LST was tested using the Mann-Kendall monotonicity test [51]. The results showed that ρ ccm between LST and all variables, except SM and AOD, showed significant monotonic increases (p < 0.05) and convergence characteristics.

Causality Detection Process between LST and Environmental Factors
To clearly express the process of causality detection between LST and environmental factors, a pixel with grassland land cover type from the PMC region was selected as an example for a detailed demonstration. To ensure more stable results, each pair of cross-mapping processes was repeated 100 times [19], and the average of the 100 results was used as the final result. Figure 3 shows the detection results between LST and each environmental factor, with the horizontal axis being the length of the time series and the vertical axis being the ρccm between LST and environmental factors.

Causality Detection Process between LST and Environmental Factors
To clearly express the process of causality detection between LST and environmenta factors, a pixel with grassland land cover type from the PMC region was selected as an example for a detailed demonstration. To ensure more stable results, each pair o cross-mapping processes was repeated 100 times [19], and the average of the 100 result was used as the final result. Figure 3 shows the detection results between LST and each environmental factor, with the horizontal axis being the length of the time series and th vertical axis being the ρccm between LST and environmental factors.  A pseudo-causal relationship between LST and environmental factors can be caused by seasonal synchronization effects. For this reason, SST tests were conducted for LST and environmental factors, except SM and AOD. The test results are shown in Figure 4, where the red circles indicate the ρ ccm between LST and environmental factors, and the boxes indicate the ρ ccm between LST and the 500 seasonal surrogate series of environmental factors. The red solid circles demonstrate that the ρ ccm between the LST and the environmental factor is significantly larger than the ρ ccm between the seasonal surrogate series and the environmental factor; that is, the causal relationship between the LST and the environmental factor is not influenced by seasonal synchronization effects. The red hollow circles indicate a pseudo-causal relationship between the LST and environmental factors caused by seasonal synchronization effects. The SST results showed a pseudo-causal relationship between PRE and LST within this pixel.
revealed that the strength of the potential causal relationships between them and LST was relatively low. The ρccm values between SM and LST, and AOD and LST decreased with increasing time series, indicating that there was no causal relationship between them and LST. The monotonic convergence of ρccm between environmental factors and LST was tested using the Mann-Kendall monotonicity test [51]. The results showed that ρccm between LST and all variables, except SM and AOD, showed significant monotonic increases (p < 0.05) and convergence characteristics.
A pseudo-causal relationship between LST and environmental factors can be caused by seasonal synchronization effects. For this reason, SST tests were conducted for LST and environmental factors, except SM and AOD. The test results are shown in Figure 4, where the red circles indicate the ρccm between LST and environmental factors, and the boxes indicate the ρccm between LST and the 500 seasonal surrogate series of environmental factors. The red solid circles demonstrate that the ρccm between the LST and the environmental factor is significantly larger than the ρccm between the seasonal surrogate series and the environmental factor; that is, the causal relationship between the LST and the environmental factor is not influenced by seasonal synchronization effects. The red hollow circles indicate a pseudo-causal relationship between the LST and environmental factors caused by seasonal synchronization effects. The SST results showed a pseudo-causal relationship between PRE and LST within this pixel. In the case of random variables, CCM and SST causality tests can be misleading because of data limitations. In this study, variables that passed the SST test were introduced into the multivariate EDM model, and this pseudo-causality was removed by examining the extent to which the variables improved the forecasting skill of the multivariate EDM model. Figure 5 shows the improvement in the forecast skill of the multivariate EDM model for LST with the addition of each environmental factor separately. The blue bars In the case of random variables, CCM and SST causality tests can be misleading because of data limitations. In this study, variables that passed the SST test were introduced into the multivariate EDM model, and this pseudo-causality was removed by examining the extent to which the variables improved the forecasting skill of the multivariate EDM model. Figure 5 shows the improvement in the forecast skill of the multivariate EDM model for LST with the addition of each environmental factor separately. The blue bars indicate the forecast skill of the multivariate EDM model for LST without the addition of environmental factors, whereas the red bars indicate the degree of improvement in the model forecast skill with the addition of environmental factors. Overall, the forecast skill of the multivariate EDM model was significantly improved by the addition of TA, NDVI, RN, ET, and WV to this pixel (p < 0.05), indicating that there was a causal relationship between the above five environmental factors and LST. RN (0.22) had the greatest improvement in the forecast skill of the model, followed by NDVI (0.17), WV (0.11), TA (0.08), and ET (0.07).
Detection of the generalized synchrony effect and the Moran effect between variables is still an indispensable step in causality detection; therefore, this study further performed the ECCM test between the LST and the five environmental factors mentioned above. The spurious driving effect of the environmental factors on LST was removed by determining the time lag length of the effect of environmental factors on LST. Figure 6 shows the results of the ECCM between the LST and each environmental factor. The time lags set in this study were four months before and after LST, and the optimal ρ ccm between the five environmental factors and LST were all negative. This demonstrated that the environmental factors had a real driving effect on LST and were not affected by the generalized synchronization effect or the Moran effect. The lag lengths also represented the lags in the driving effect of environmental factors on LST.
indicate the forecast skill of the multivariate EDM model for LST without the addition of environmental factors, whereas the red bars indicate the degree of improvement in the model forecast skill with the addition of environmental factors. Overall, the forecast skill of the multivariate EDM model was significantly improved by the addition of TA, NDVI, RN, ET, and WV to this pixel (p < 0.05), indicating that there was a causal relationship between the above five environmental factors and LST. RN (0.22) had the greatest improvement in the forecast skill of the model, followed by NDVI (0.17), WV (0.11), TA (0.08), and ET (0.07). Detection of the generalized synchrony effect and the Moran effect between variables is still an indispensable step in causality detection; therefore, this study further performed the ECCM test between the LST and the five environmental factors mentioned above. The spurious driving effect of the environmental factors on LST was removed by determining the time lag length of the effect of environmental factors on LST. Figure 6 shows the results of the ECCM between the LST and each environmental factor. The time lags set in this study were four months before and after LST, and the optimal ρccm between the five environmental factors and LST were all negative. This demonstrated that the environmental factors had a real driving effect on LST and were not affected by the generalized synchronization effect or the Moran effect. The lag lengths also represented the lags in the driving effect of environmental factors on LST.
In summary, after performing CCM causality tests between LST and environmental factors, and conducting pseudo-causality tests, such as SST, multivariate EDM, and ECCM tests, we identified TA, NDVI, RN, ET, and WV as drivers of LST change within this pixel.
Once the drivers of LST were identified, the direction and strength of the effects of the drivers on LST were analyzed using multivariate scenario analysis. Figure 7 shows the degree of response of the LST to small changes in the five drivers within this pixel, with the horizontal axis showing the drivers and the vertical axis indicating the degree of response of the LST. The results revealed that TA, NDVI, and RN had a positive driving effect on LST; that is, LST increased as the driving factor increased. Moreover, ET and WV had a negative driving effect on LST, that is, LST decreased as the driving factor increased. In addition, RN had the strongest effect on LST with a mean value of 0.58, followed by NDVI with a mean value of 0.13. TA had a relatively weak positive driving effect on LST with a mean value of 0.10. ET and WV had weak negative effects on LST with mean values of −0.04 and −0.02, respectively.  In summary, after performing CCM causality tests between LST and environmental factors, and conducting pseudo-causality tests, such as SST, multivariate EDM, and ECCM tests, we identified TA, NDVI, RN, ET, and WV as drivers of LST change within this pixel.
Once the drivers of LST were identified, the direction and strength of the effects of the drivers on LST were analyzed using multivariate scenario analysis. Figure 7 shows the degree of response of the LST to small changes in the five drivers within this pixel, with the horizontal axis showing the drivers and the vertical axis indicating the degree of response of the LST. The results revealed that TA, NDVI, and RN had a positive driving effect on LST; that is, LST increased as the driving factor increased. Moreover, ET and WV had a negative driving effect on LST, that is, LST decreased as the driving factor increased. In addition, RN had the strongest effect on LST with a mean value of 0.58, followed by NDVI with a mean value of 0.13. TA had a relatively weak positive driving effect on LST with a mean value of 0.10. ET and WV had weak negative effects on LST with mean values of −0.04 and −0.02, respectively.

Result of CCM and SST Tests of China
Causality tests were conducted on a pixel-by-pixel basis to examine the causal relationships between LST and each environmental factor across China. The spatial distribution of ρ ccm between the LST and each environmental factor, as shown in Figure 8, was the result of the CCM and SST tests. The right-hand side of each panel corresponds to the latitudinal distribution of the mean value of ρ ccm . Overall, the ρ ccm of TA and RN were above 0.8 and 0.7, respectively, followed by NDVI and ET. However, the spatial variability between them was large, with variations ranging from 0.1 to 0.98 and 0.2 to 0.95, respectively. The ρ ccm of SM, PRE, AOD, and WV were relatively small.
From the perspective of latitudinal distribution, the ρ ccm of TA increased with increasing latitude and was less than the mean value of 0.94 south of 37 • N, and greater than the mean value of 0.94 between 37 • N and 53 • N. It saw a decline north of 53 • N. The ρ ccm of NDVI also increased with increasing latitude, but increased rapidly between 30 • N and 35 • N, and then decreased between 35 • N and 40 • N. This may be because the region between 30-35 • N is the highest altitude Qinghai-Tibet Plateau region, while 35 • -40 • N falls within the largest desert in China. The ρ ccm of the SM varied less with latitude, fluctuating around a mean value of 0.43. The ρ ccm of RN had similar distribution characteristics to TA with respect to latitude, with a mean value of 0.91. The ρ ccm of the PRE was generally lower south of 32 • N and north of 53 • N, as compared to the other latitudinal zones. The ρ ccm of AOD was significantly higher than the mean between 38 • N and 43 • N. This is because the main surface cover type in this region was bare ground with little rain, high dust content in the air, and increased aerosol impact on LST. The ρ ccm of ET gradually increased from low to high latitudes, with a mean value of 0.70. The ρ ccm value of WV also increased with increasing latitude, but with a smaller mean value of 0.31.
Causality tests were conducted on a pixel-by-pixel basis to examine the causal rela tionships between LST and each environmental factor across China. The spatial distribu tion of ρccm between the LST and each environmental factor, as shown in Figure 8, was th result of the CCM and SST tests. The right-hand side of each panel corresponds to th latitudinal distribution of the mean value of ρccm. Overall, the ρccm of TA and RN wer above 0.8 and 0.7, respectively, followed by NDVI and ET. However, the spatial varia bility between them was large, with variations ranging from 0.1 to 0.98 and 0.2 to 0.95 respectively. The ρccm of SM, PRE, AOD, and WV were relatively small. From the perspective of latitudinal distribution, the ρccm of TA increased with in creasing latitude and was less than the mean value of 0.94 south of 37°N, and greater than the mean value of 0.94 between 37°N and 53°N. It saw a decline north of 53°N. The ρccm o NDVI also increased with increasing latitude, but increased rapidly between 30°N and 35°N, and then decreased between 35°N and 40°N. This may be because the region be tween 30-35°N is the highest altitude Qinghai-Tibet Plateau region, while 35°-40°N fall within the largest desert in China. The ρccm of the SM varied less with latitude, fluctuating around a mean value of 0.43. The ρccm of RN had similar distribution characteristics to TA with respect to latitude, with a mean value of 0.91. The ρccm of the PRE was generally  Figure 9 presents the statistics of ρ ccm between LST and each environmental factor across climate zones. The results showed that the ρ ccm values of TA and RN were still relatively strong compared to other factors in the same climate zone. TA and RN had the strongest ρ ccm in TCC regions with mean values of 0.98 and 0.97, respectively, while they had the weakest ρ ccm in TRMC with mean values of 0.89 and 0.90, respectively. For NDVI and ET, the strongest mean values of ρ ccm were found in the TMC region, with mean values of 0.84 and 0.80, respectively. This may be due to the significant seasonal variation in vegetation in the TMC region, which strengthens the influence of NDVI on LST. In addition, the warmer and wetter monsoon climate enhanced the role of ET. However, the ρ ccm values of NDVI and ET were the weakest in the TRMC region, with mean values of 0.38 and 0.54, respectively. This suggested that the strength of the potential causal relationships between NDVI and LST as well as ET and LST were limited by hydrothermal conditions. ues of 0.84 and 0.80, respectively. This may be due to the significant seasonal variation in vegetation in the TMC region, which strengthens the influence of NDVI on LST. In addi tion, the warmer and wetter monsoon climate enhanced the role of ET. However, the ρccm values of NDVI and ET were the weakest in the TRMC region, with mean values of 0.38 and 0.54, respectively. This suggested that the strength of the potential causal relationships between NDVI and LST as well as ET and LST were limited by hydrothermal conditions.   Figure 10 demonstrates the results of the multivariate EDM test between the environmental factors and LST after detection using CCM and SST. It can be seen that the results of multivariate EDM are presented in two categories: the forecast skill of the model with improvement (blue) and without improvement (yellow). For regions where the forecast skill of the multivariate EDM was not improved, causal relationship between the environmental factor and LST was considered to be pseudo-causal.

Pseudo-Causality Detection by Multivariate EDM
As can be seen from Figure 10, the forecast skill of the multivariate EDM with the addition of RN was significantly improved over the full domain, with only a sporadic area in the southern part of the Tibetan Plateau showing no improvement, accounting for less than 1%. Apart from RN, the forecast skill of the model after the addition of TA, PRE, ET, or WV also showed remarkable improvement; the regions that did not show betterment accounted for less than 10%. Pseudo-causality between NDVI and LST accounted for 13.29% of the total, mainly in the northern high latitudes where vegetation was sparse. Multivariate EDM detected a high proportion of pseudo-causal relationships between AOD and LST, at 25.26%, which was spread over most of the eastern region. SM had the highest pseudo-causal relationship with LST (43.33%), occupying almost all of northern China. Figure 10 demonstrates the results of the multivariate EDM test between the ronmental factors and LST after detection using CCM and SST. It can be seen tha results of multivariate EDM are presented in two categories: the forecast skill o model with improvement (blue) and without improvement (yellow). For regions w the forecast skill of the multivariate EDM was not improved, causal relationship tween the environmental factor and LST was considered to be pseudo-causal. As can be seen from Figure 10, the forecast skill of the multivariate EDM with the tion of RN was significantly improved over the full domain, with only a sporadic area i southern part of the Tibetan Plateau showing no improvement, accounting for less tha Apart from RN, the forecast skill of the model after the addition of TA, PRE, ET, or WV showed remarkable improvement; the regions that did not show betterment accounte less than 10%. Pseudo-causality between NDVI and LST accounted for 13.29% of the mainly in the northern high latitudes where vegetation was sparse. Multivariate EDM tected a high proportion of pseudo-causal relationships between AOD and LST, at 25 which was spread over most of the eastern region. SM had the highest pseudo-causal tionship with LST (43.33%), occupying almost all of northern China.  Figure 11 provides an overview of the results of ECCM detection between LST and environmental factors. Numerous positive time lags were detected for all factors due to generalized synchrony and Moran effects. Of these, the positive lag between TA and LST accounted for just 10.16%. Nearly 50% of these were distributed in the PMC region, probably because the PMC region was at an extremely high altitude and had the lowest air temperature in the study area. Therefore, its causal relationship with LST was not significant. The positive lag between RN and LST was 24.9%; approximately 20% of the positive lags were distributed in the TCC and TMC regions. With the exception of TA and RN, the positive lag between LST and other environmental factors was evenly distributed across different climatic zones, accounting for over 30%. In particular, two factors, AOD and SM, accounted for approximately 50% of the positive lags. Thus, the causality of PRE and AOD with LST was not only weak, but at least half of them were false. However, a positive time lag detected between environmental factors and LST only indicates that environmental factors have no driving effect on LST. They do not negate the absence of a driving effect of LST on environmental factors.

Pseudo-Causality Detection by ECCM
probably because the PMC region was at an extremely high altitude and had th air temperature in the study area. Therefore, its causal relationship with LST significant. The positive lag between RN and LST was 24.9%; approximately 20 positive lags were distributed in the TCC and TMC regions. With the exceptio and RN, the positive lag between LST and other environmental factors was eve tributed across different climatic zones, accounting for over 30%. In particular, tors, AOD and SM, accounted for approximately 50% of the positive lags. Thus, sality of PRE and AOD with LST was not only weak, but at least half of them we However, a positive time lag detected between environmental factors and LST dicates that environmental factors have no driving effect on LST. They do not ne absence of a driving effect of LST on environmental factors. After pseudo-causal elimination with SST, multivariate EDM, and ECCM, portion of each environmental factor as a driver of the LST was determined. A in Figure 12, the highest percentage of driving factors among the environmenta After pseudo-causal elimination with SST, multivariate EDM, and ECCM, the proportion of each environmental factor as a driver of the LST was determined. As shown in Figure 12, the highest percentage of driving factors among the environmental factors at the national scale was TA, accounting for 50.11%. ρ ccm had a mean value of 0.96. From different climate zones, the mean value of ρ ccm between TA and LST, as well as the percentage as a driving factor, showed a gradual decrease from north to south, indicating that LST was influenced by TA in its scope and intensity, both of which gradually weakened from north to south. The ρ ccm of TA within the TCC was the largest, with a mean value of 0.98, and the percentage of driving factor was 70%. The mean value of ρ ccm within the TRMC was the smallest (0.90), with 11.08% as the driving factor. In addition to TA, RN also had a relatively high share of 46.9%, with a mean ρ ccm of 0.92. The climatic zone with the largest ρ ccm of RN was TCC with a mean value of 0.96 and a driver share of 46.91%, while the climatic zone with the largest driver share was PMC (57.74%), with a mean ρ ccm of 0.92. The proportion of NDVI as a driver of LST at the national scale was 34.28%, with a mean ρ ccm of 0.65. TMC was the climate zone with the largest ρ ccm of NDVI, with a mean value of 0.84 and a percentage of drivers of 45.87%. ET as a driver of LST was relatively low at 25.4%, but ρ ccm was higher than NDVI with a mean value of 0.74. The climate zone with the highest proportion of ET as a driver of LST was TMC, with a proportion of 40.54% and a mean value of 0.81 for ρ ccm . The percentages of drivers in SM, AOD, PRE, and WV were 16.43%, 19.71%, 30.14%, and 30.36%, respectively, and with the mean ρ ccm of 0.48, 0.37, 0.36, and 0.34, respectively. at the national scale was TA, accounting for 50.11%. ρccm had a mean value of 0.96. From different climate zones, the mean value of ρccm between TA and LST, as well as the per centage as a driving factor, showed a gradual decrease from north to south, indicating that LST was influenced by TA in its scope and intensity, both of which gradually weakened from north to south. The ρccm of TA within the TCC was the largest, with mean value of 0.98, and the percentage of driving factor was 70%. The mean value of ρccm within the TRMC was the smallest (0.90), with 11.08% as the driving factor. In addition to TA, RN also had a relatively high share of 46.9%, with a mean ρccm of 0.92. The climat ic zone with the largest ρccm of RN was TCC with a mean value of 0.96 and a driver shar of 46.91%, while the climatic zone with the largest driver share was PMC (57.74%), with a mean ρccm of 0.92. The proportion of NDVI as a driver of LST at the national scale wa 34.28%, with a mean ρccm of 0.65. TMC was the climate zone with the largest ρccm o NDVI, with a mean value of 0.84 and a percentage of drivers of 45.87%. ET as a driver o LST was relatively low at 25.4%, but ρccm was higher than NDVI with a mean value o 0.74. The climate zone with the highest proportion of ET as a driver of LST was TMC with a proportion of 40.54% and a mean value of 0.81 for ρccm. The percentages of driver in SM, AOD, PRE, and WV were 16.43%, 19.71%, 30.14%, and 30.36%, respectively, and with the mean ρccm of 0.48, 0.37, 0.36, and 0.34, respectively.

Effect of Driving Factors on LST
Having identified the drivers of LST over the study period, we analyzed the driving effect of the drivers. Figure 13 shows the spatial distribution of the results of the scenario analysis between each driver and LST. The red areas indicate a positive effect of the drivers on LST and the blue areas indicate a negative effect of the drivers on LST. The results revealed that there was a positive effect on LST at the national level, as shown in Figure 13a,d. A similar relationship between TA, RN, and LST was also reported in previous studies [52][53][54]. The positive effects of NDVI on LST were mostly found in the TCC and PMC regions ( Figure 13b); however, both climate regions had restricted energy access. The negative effect of NDVI on LST was mostly observed in the TMC, sub_TRMC, and TRMC regions. This result reaffirmed the conclusion of Karnieli et al. [55] that the interaction between LST and NDVI is positive when energy is limited (high latitude or high altitude) and negative when energy is available [56]. In addition, the positive effect of ET on LST was mainly concentrated in the eastern part of the TCC region, most of the PMC region, and central part of the sub_TRMC region. Sun et al. [2] stated that the limits of TA and SM should be considered when analyzing the relationship between ET and LST. However, in this study, the interaction between ET and LST were related to land cover type. AOD had a predominantly positive effect on LST, whereas the regions that showed a negative effect were mainly located in the western part of the TCC region. Studies have shown that the impact of aerosols on climate depends not only on the quantity of aerosols but also on their optical properties. Aerosols affect LST by influencing the solar radiation reaching the surface. Scattering aerosols having a cooling effect as they reduce the energy received at the surface, whereas absorbing aerosols have a warming effect [55]. In different climatic regions, there were significant differences in the manner and intensity of the effect of the driver (Figure 14). At a national level, the positive effect of TA on LST was 46.42%, whereas the negative effect was only 0.69%. In terms of the mean value of the intensity of the effect TA on LST, it was maximum in the TMC region with a mean value of 0.97, in minimum in the TRMC region with a mean value of 0.35. This indicated that the effect of TA on LST was stronger at higher latitudes. The strength of RN's effect on LST was second only to that of TA, with 41% positive and 1.4% negative effects. Regarding the mean value of the effect of RN on LST, the climatic region with the strongest effect of RN on LST was TMC with a mean value of 0.69; the weakest was TRMC with a mean value of 0.22. The positive effect of NDVI on LST was 22.27% and the negative effect was 8.56%. Of these, TCC, TMC, PMC, and sub_TRMC were positive and the highest within PMC, with a mean value of 0.15, while in TRMC, they were mainly negative with a mean value of −0.01. On LST, ET had a positive effect of 15.11% and a negative effect of 9.93%, with a greater intensity within PMC and sub_TRMC, with mean values of 0.19 and 0.17, respectively. In different climatic regions, there were significant differences in the manner and intensity of the effect of the driver (Figure 14). At a national level, the positive effect of TA on LST was 46.42%, whereas the negative effect was only 0.69%. In terms of the mean value of the intensity of the effect TA on LST, it was maximum in the TMC region with a mean value of 0.97, in minimum in the TRMC region with a mean value of 0.35. This indicated that the effect of TA on LST was stronger at higher latitudes. The strength of RN's effect on LST was second only to that of TA, with 41% positive and 1.4% negative effects. Regarding the mean value of the effect of RN on LST, the climatic region with the strongest effect of RN on LST was TMC with a mean value of 0.69; the weakest was TRMC with a mean value of 0.22. The positive effect of NDVI on LST was 22.27% and the negative effect was 8.56%. Of these, TCC, TMC, PMC, and sub_TRMC were positive and the highest within PMC, with a mean value of 0.15, while in TRMC, they were mainly negative with a mean value of −0.01. On LST, ET had a positive effect of 15.11%

Type of Causal Relationship between LST and Drivers
Recent studies have focused on the relationship between single variables and LST [10,11,14,57]. The identification of the driving role of multiple variables on LST and their dominant factors has rarely been studied. Furthermore, correlation analysis has been used in the methods used to study the relationship between LST and influencing factors, but it needs to be emphasized that correlation does not mean that a causal relationship exists. This study attempted to explore the dominant drivers of LST in terms of the causal relationship between multiple factors and LST. The CCM method, which is suitable for causality detection in nonlinear dynamical systems, was used to detect the causal relationship between each environmental factor and LST. The CCM test results were verified by SST, multivariate EDM, and ECCM, which showed that TA and RN were strong drivers of LST and had a much higher positive effect than other factors. Furthermore, NDVI and ET had relatively strong driving effects on LST, but they had both positive and negative effects in different regions.
In addition, the CCM method was used to identify whether there was a causal relationship between environmental factors and LST, but it could not determine whether the driving effect was unidirectional or bidirectional. To determine the interaction pattern between each driver and LST, this study used the CCM method to test the driving effect of LST on environmental factors and performed a pseudo-causality test. The interactions between LST and environmental factors were classified into three categories: environmental factors driving LST, LST driving environmental factors, and mutual drive. Figure 15 show the spatial distribution of the pattern of causal interaction between LST and drivers and the share of each type, respectively. Total of the three causality types between RN and LST accounted for 79.38% of the country and 20.62% had no causality (Figure 15d). Among them, bidirectional causality accounted for 48.49% and was widely distributed, while the type of LST-driven RN accounted for 20.23% and was mainly distributed in the TCC and southern part of TMC, where RN was lower. Second to RN, the areas with causality between TA and LST accounted for 78.04% of the whole country land area, while 21.96% of the areas had no causality (Figure 15a). The type of causality was mainly bidirectional, accounting for 63.36% of the country, and was mainly distributed in the TCC, TMC, and PMC regions. The causality between NDVI and LST was 34.18%, 16.31%, and 16.34% for bidirectional, LST-driven NDVI, and NDVI-driven LST, respectively, and 33.18% had no causality (Figure 15b). PRE was the only variable among the selected environmental factors that was dominated by reverse causality, that is, LST-driven PRE type causality was dominant (31.32%), and mainly in the TCC and PMC areas, which was higher than the proportion of bidirectional causality at 24.79% (Figure 15e). In addition, the WV-driven LST causality type accounted for the largest proportion of the three causality types (25.21%). For SM and AOD, no causality has a large percentage of the country, with 68.81% and 57.9%, respectively.
Overall, the causality between LST and most drivers was predominantly bidirectional. On the other hand, few drivers, such as PRE and WV, were predominantly unidirectional. Overall, the causality between LST and most drivers was predominantly bidirect On the other hand, few drivers, such as PRE and WV, were predominantly unidirection

Bidirectional Causality between LST and Drivers
For the bidirectional causality type, there is a need to understand the strong rection. To this end, we calculated the difference (Δρccm) between ρccm in the positiv rection (factor driving LST) and ρccm in the negative direction (LST driving factor).

Bidirectional Causality between LST and Drivers
For the bidirectional causality type, there is a need to understand the stronger direction. To this end, we calculated the difference (∆ρ ccm ) between ρ ccm in the positive direction (factor driving LST) and ρ ccm in the negative direction (LST driving factor). Figure 16 shows the variation in ∆ρ ccm with latitude. Overall, all variables were dominated by negative differences, except for TA, which was dominated by positive differences. In the interaction between TA and LST, the intensity of the TA driving LST north of 37 • N was higher than that of the LST driving TA. The difference between 37 • N and Remote Sens. 2022, 14, 3280 22 of 26 33 • N was close to 0, indicating that the strength of the interaction between LST and TA was similar in this interval. Negative peaks occurred between 33 • N and 27 • N, and fluctuations in the difference were greater south of 27 • N, which may be due to the small amount of data in this region. Generally, the difference in the strength of the causality between TA and LST in both directions was small, with ∆ρ ccm ranging from −0.05 to 0.1. The causality between NDVI and LST was dominated by the LST-driven NDVI direction throughout the country. The absolute value of ∆ρ ccm decreased gradually from south to north, indicating that LST drove NDVI more in the south. The difference in the strength of the two-way causality between RN and LST was also relatively small, dominated by LST-driven RN. There was a trend of stronger RN driving LST south of 25 • N, with a negative peak between 25 • N and 40 • N. This decreased rapidly because of the extremely high altitude of the Tibetan Plateau in this region. The ∆ρ ccm between AOD and LST only showed a positive peak between 30 • N and 35 • N, with negative values in all other latitudinal zones, as the Beijing-Tianjin-Hebei region with high AOD values was distributed in this latitudinal zone. The strength of AOD driving LST was higher than that of LST driving AOD. The relationship between ET and LST was also dominated by LST driving ET, and ∆ρ ccm was relatively stable at all latitudes.  Figure 16 shows the variation in Δρccm with latitude. Overall, all variables were dominated by negative differences, except for TA, which was dominated by positive differences. In the interaction between TA and LST, the intensity of the TA driving LST north of 37°N was higher than that of the LST driving TA. The difference between 37°N and 33°N was close to 0, indicating that the strength of the interaction between LST and TA was similar in this interval. Negative peaks occurred between 33°N and 27°N, and fluctuations in the difference were greater south of 27°N, which may be due to the small amount of data in this region. Generally, the difference in the strength of the causality between TA and LST in both directions was small, with Δρccm ranging from −0.05 to 0.1. The causality between NDVI and LST was dominated by the LST-driven NDVI direction throughout the country. The absolute value of Δρccm decreased gradually from south to north, indicating that LST drove NDVI more in the south. The difference in the strength of the two-way causality between RN and LST was also relatively small, dominated by LST-driven RN. There was a trend of stronger RN driving LST south of 25°N, with a negative peak between 25°N and 40°N. This decreased rapidly because of the extremely high altitude of the Tibetan Plateau in this region. The Δρccm between AOD and LST only showed a positive peak between 30°N and 35°N, with negative values in all other latitudinal zones, as the Beijing-Tianjin-Hebei region with high AOD values was distributed in this latitudinal zone. The strength of AOD driving LST was higher than that of LST driving AOD. The relationship between ET and LST was also dominated by LST driving ET, and Δρccm was relatively stable at all latitudes.  In summary, even though there was a bidirectional causal relationship between LST and each environmental factor, the dominant causal relationship was still in the direction of LST driving environmental factors.
In this study, CCM, SST, multivariate EDM, and ECCM were used to identify the true drivers of LST among environmental factors and to explore the interaction patterns between LST and drivers through multivariate scenario analysis. However, the complexity of surface eco-physical systems prevented a complete list of all variables associated with LST, and the diversity of data sources introduced inevitable errors in the experimental results. Therefore, a complete understanding of LST can be obtained from an integrated treatment of these constraints, which will provide a basis for future research.

Conclusions
This study explored the driving role of environmental factors (TA, NDVI, SM, RN, PRE, AOD, ET, and WV) and LST during 2003-2018 in China using the CCM, SST, multivariate EDM, and ECCM methods.
All environmental factors had driving effects on LST. However, after the removal of pseudo-causality and the comparison of the driving strength of each variable, TA and RN were concluded to be the dominant factors with high driving strengths, followed by NDVI and ET.
All drivers had positive as well as negative effects on LST, except for TA and RN, which only had positive effects on LST. The intensity of the effects of TA and RN on LST showed signs of weakening from high northern latitudes to low latitudes, indicating that LST was more strongly driven by TA and RN in regions with lower energy. The positive effects of NDVI and ET on LST were mainly distributed in the TCC and PMC zones, while the negative effects were mainly distributed in the TMC, sub_TRMC, and TRMC regions. This indicated that the direction of NDVI and ET effects on LST were limited by hydrothermal conditions. Overall, the positive effects were more significant at high latitudes and altitudes.
There are three types of causal relationships between LST and each environmental factor, that is, factor-driven LST, LST-driven factor, and bidirectional drive. However, the bidirectional causal relationship between LST and the driver is dominant. In addition, the categorization of causal relationships between LST and factors revealed that LST driving factors were dominant in quantity, except for SM and WV. In terms of driving strength, the driving effect of LST on factors was stronger than the driving strength of factors on LST, except for TA.
Our results quantify the driving effects of environmental factors on LST; however, future studies should explore the relationships between surface type, elevation, and slope orientation to fully understand the dynamics of LST.