*3.2. LST, Vegetation, and Albedo Trends during 2000–2021 at Annual and Seasonal Scales* 3.2.1. LST Trends at Different Timescales

The spatial pattern of LST trends recorded at different timescales was more concentrated than that for the NDVI data. Alongside being common in winter, the LST warming trend on the TP from 2000 to 2021 was most prominent in the southern TP (Figure 7a–e). The LST warming rate was significantly higher in summer and autumn than during winter and spring, and it was mainly concentrated in the southwestern TP (Figure 7c,d). The LST cooling trend was most prominent on the northern TP on an annual scale and during spring, summer, and autumn (Figure 7a–d). During the winter, cooling mostly occurred in the southwestern TP. The spatial distribution of significance trends (*p* < 0.05) showed that the most significant decrease in LST occurs at annual scale and in spring (Figure 7f,j), with 7.06% and 7.93% of all pixels recording these changes, respectively. Significant LST increases (*p* < 0.05) in autumn were represented by 10.60% of all pixels (Figure 7i), which covered a much larger area than that for spring (5.01%; Figure 7g), summer (6.71%; Figure 7h), and winter (1.42%; Figure 7j).

**Figure 7.** MODIS-based LST trends during 2000–2021 at annual and seasonal scales. (**a**–**e**) denote the LST trends for 2000–2021 in annual, spring, summer, autumn, and winter, respectively; (**f**–**j**) denote the LST trends with *p* < 0.05 for 2000–2021 in annual, spring, summer, autumn, and winter, respectively.

### 3.2.2. NDVI Trends at Different Timescales

Spatial variations in the surface parameters across the study region revealed that MODIS-based NDVI data show significantly more greening trends than GIMMS-based NDVI data (Figure 8). These NDVI data showed that the eastern TP records the largest greening area during 2000–2021 at all timescales (Figure 8a–e), and the greening rate in the spring is significantly higher than for the other timescales (Figure 8b). The southern TP showed the highest NDVI browning rate (Figure 8a–d), which was highest in autumn out of all the considered timescales (Figure 8d). The NDVI data recorded an increasing trend on an annual scale, with 40.16% of pixels showing a significant increase and only 3.56% showing a significant decrease trend (Figure 8f, *p* < 0.05). Approximately 44.62% of the NDVI pixels showed a significant increase in the spring, and 2.53% significantly decreased (Figure 8g). Furthermore, 38.84% of the NDVI pixels showed significant greening (*p* < 0.05), and 3.49% showed a considerable browning during the summer (Figure 8h). An approximately equal proportion of NDVI pixels showed browning (21.13%) and greening (21.35%) in autumn (Figure 8i). Finally, 52.56% and 2.24% of the NDVI pixels recorded greening and browning in winter, respectively, with these values being significantly higher than their equivalents during the other seasons (Figure 8j).

### 3.2.3. Albedo Trends at Different Timescales

Albedo trends at all timescales showed relatively few spatial hotspots except for in winter (Figure 9a–e). The fastest albedo growth rates during 2000–2021 occurred in spring and winter (Figure 9b,e); however, the albedo trends in summer and autumn were less pronounced (Figure 9c,d). The decreasing rate of albedo was similar to its increasing rate in the temporal pattern. A significant increase in albedo was observed in the western plateau region at all timescales, and a significant decrease is found in the south and northeast plateau (Figure 9f–j, *p* < 0.05). In terms of areas distribution, a relatively high percentage of pixels showed a significant increase during spring (7.78%; Figure 9g) and winter (4.92%; Figure 9j), and they were notably higher than the values of 0.49% in summer and 0.38% in autumn. This effect may have been due to snow accumulation in spring and winter. The largest areas with significant decreasing trends in albedo were monitored in summer and autumn (Figure 9h,i, *p* < 0.05). These areas comprised 26.30% and 16.36% of the total number of pixels in the study region, respectively, and these changes might potentially be linked to increased vegetation coverage during the growing season.

### 3.2.4. Vegetation Impacts on LST at Different Time Scales

A statistical method was used to quantify the individual effect of vegetation on LST, and we simultaneously considered the varying albedo (see Section 3.2.3) due to its significant role in the vegetation–LST relationship [46,47]. Specifically, we employed the detrending method and partial correlation analysis to analyze the individual effects of the NDVI and albedo on LST. First, we examined the linear trends of the NDVI and albedo (see the Supporting Information S3 for details). The results suggested that the NDVI of the TP shows significant temporal trends at different timescales, and the summer albedo also shows significant temporal trends (Table 2). Consequently, we filtered the temporal trends of the indicated six variables using the detrending method. Second, we fixed NDVI (albedo) to analyze the correlation between albedo (NDVI) and LST. The results showed that the partial correlation coefficients of the NDVI and albedo with LST are significant except in summer (Table 3). This may be due to the saturation effect of the summer vegetation and albedo contributions on the LST [48], resulting in their insignificant trend of contribution. The individual contribution of albedo to LST was larger than that of NDVI at all timescales, especially in spring and autumn. This may be attributed to the fact that (1) the individual contribution of albedo to LST includes both the indirect contribution of vegetation altering albedo [49] and the direct contribution of albedo; and (2) an advanced vegetation growing season and a delayed vegetation ending season can alter surface albedo strongly [50,51], causing the contribution of surface albedo to be larger in spring and autumn. In terms

of winter, faded vegetation and deepened snow cover would reduce the contribution of vegetation but increase the contribution of albedo [47].

**Figure 8.** MODIS-based NDVI trends during 2000–2021 at annual and seasonal scales. (**a**–**e**) denote the NDVI trends for 2000–2021 in annual, spring, summer, autumn, and winter, respectively; (**f**–**j**) denote the NDVI trends with *p* < 0.05 for 2000–2021 in annual, spring, summer, autumn, and winter, respectively.

**Figure 9.** MODIS-based albedo trends during 2000–2021 at annual and seasonal scales. (**a**–**e**) denote the albedo trends for 2000–2021 in annual, spring, summer, autumn, and winter, respectively; (**f**–**j**) denote the albedo trends with *p* < 0.05 for 2000–2021 in annual, spring, summer, autumn, and winter, respectively.


**Table 2.** Significance of temporal trends in NDVI and albedo over 2000–2021.

**Table 3.** Partial correlation coefficients and significance of NDVI (albedo) with LST after fixing for albedo (NDVI).


Note that we utilized a statistical method to elucidate the relationship between vegetation, albedo, and LST, but we acknowledge its limitations due to the absence of a detailed surface radiation balance analysis. Achieving an accurate decomposition of surface radiation components would require further discussion beyond the scope and word limit of this paper. Therefore, for the purpose of this study, we considered the statistical method as a reasonable proxy for interpreting the individual impacts of vegetation and albedo on LST at various timescales. The individual impacts of vegetation and albedo to LST using surface radiative balance methods is expected in the near future. In addition, vegetation data are considered in producing LST products, yet it is difficult to estimate the contribution of vegetation data due to a lack of robust methods. Therefore, the vegetation contribution to LST estimation may have some uncertainties in the dense vegetation growth of the southeastern TP.

### **4. Discussion**

Vegetation–temperature relationships can be affected by estimation methods, data sources, and land cover types. We therefore discussed the breakpoint method for vegetation trends, the differences between GIMMS and MODIS NDVI datasets, and the impacts of land cover types on annual NDVI trends.

Firstly, we analyzed vegetation–temperature relationships during different warming phases. In terms of methodology, the relationship between vegetation and climate variables on the TP presented remarkable phase changes. However, the phase division of such a relationship was usually subject to the limitation of the time span of remote sensing observed vegetation datasets and the study period. For example, a 5-year time step was used to investigate the relationship between climatic and non-climatic factors and vegetation on the TP during 1980–2010 [40], and a visual graphical approach was used to determine the breakpoints in TP vegetation during 1982–2002 [41] as well as artificial time segmentation to investigate the TP vegetation feedback to climatic factors [39,46]. These methods may lead to bias in the vegetation response analysis. We therefore used the breakpoints in temperature change as the reference time points when analyzing the TP vegetation changes at different phases. We found that the breakpoints of TP temperature during 1982–2013 at different time scales largely occurred between 1994 and 1998 (Figures 2–6). This result is not only closer to direct NDVI segmentation investigation on the TP [52], but it is also closer to the artificial breakpoint of 2000 years in some investigations of TP NDVI response analysis [39,46]. Nevertheless, our results strongly imply that the artificial abrupt points later than 2000 may bring greater bias in TP NDVI response analysis. Our study also found that a common warming trend occurred at all timescales during 1982–2013, and the breakpoints of the TP were earlier in spring and autumn than those annually and in the summer and winter (Figures 2–6). This indicates that the growth and end season of vegetation may advance the warming breakpoint of the TP [50]. With the TP warming, the second phase of NDVI increase was significant at all timescales except winter, implying potential temperature accumulation induced by the first warming phase contributing to the vegetation growth in the second phase.

Secondly, our investigation demonstrated a significant NDVI trend difference between the GIMMS and MODIS datasets. Overall, NDVI trends were greater in MODIS data than in GIMMS, and this finding was consistent with [46,53]. The GIMMS NDVI data did not show a significant increase in 2000–2013, which was also found in [46,54]. This may be attributed to the wider NIR band [55] and the imperfect atmospheric effects processing technique for GIMMS data [56]. Fortunately, MODIS data are usually more reliable than GIMMS data and fill the GIMMS data gap well [57]. For example, our results indicated a significant increase in MODIS NDVI trends during 2000–2021 (Figure 8). Nevertheless, both GIMMS and MODIS products showed a common trend of vegetation changes on the TP (Figures 3a and 8a). However, given the advantages of the long time series of GIMMS data and the moderate accuracy of MODIS data, it is of great significance to integrate multiple datasets to analyze the future vegetation and climatic factors on the TP.

Thirdly, given the recent warm and wet trends on the TP [31,58], we further estimated the relative impact of NDVI changes in different land cover types on annual NDVI trend on the TP. Based on the MODIS land cover product, we detected annual land cover data in 2001 and 2021 and produced change pixels and unchanged land cover types on the TP. The results showed that 89% of the land cover types on the TP have remained unchanged over the 2001–2021 period (Figure 10). Given that the unchanged pixels cover the majority of the TP, a stepwise backward method based on these areas was employed to filter the intended land cover types and then evaluate their relative impact on the estimation of NDVI trends of TP. The results suggested that only bare land has a relative impact of 9.746% on the TP NDVI trend, but the relative impact of other land cover types on the TP NDVI trend does not exceed 5% (Table 4). This implied that the NDVI growth of bare land contributes significantly to the TP vegetation. In addition, we also found that the mean NDVI value of the TP reaches the highest after filtering the bare land compared to other cover types (Figure 11). Accordingly, we suggest that future NDVI trend and intensity studies on the TP should pay more attention to vegetation changes in bare land.

**Table 4.** Relative impact of land cover on annual NDVI trends on the Tibetan Plateau.




Notes: Null means no filtered land cover type.

**Figure 10.** Land cover change detection on the Tibetan Plateau during 2001–2021.

**Figure 11.** Mean NDVI distribution for 2000–2021 after removing different land cover types. L1, L2, L3, L4, L5, L6, L7, L8, L9, L10, L11, L12, L13, L14, L15, L16, and L17 represent evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, mixed forests, closed shrublands, open shrublands, woody savannas, savannas, grasslands, permanent wetlands, croplands, urban and built-up lands, cropland/natural vegetation, mosaics, permanent snow and ice, barren, and water bodies, respectively. Null means no filtered land cover type.

### **5. Conclusions**

This study aimed to investigate the differences in vegetation changes during various warming phases on the Tibetan Plateau (TP) as well as the trends and patterns of NDVI at different timescales and its influence on LST over the past two decades. First, the warming trend on the TP occurred at different timescales, corresponding to breakpoints that took place in 1994–1998, with earlier occurrences observed in spring and autumn. Secondly, as the warming of the TP continued, the vegetation in the second warming phase exhibited a larger greening area. These findings suggested that the time of temperature breakpoints may have an impact on the area of vegetation increase, and NDVI changes are more sensitive to air temperature during the recent warming phase. Third, MODIS data highlighted more greening compared to GIMMS data. Using MODIS data, we also found the fastest NDVI increase trend in spring and the fastest LST warming in autumn. The partial correlation analysis indicated that NDVI has a significant negative impact on LST during the annual scale and autumn while also having a significant positive impact on LST during spring. This suggests that the contribution of NDVI to LST varies across different timescales since 2000. To summarize, our findings systematically uncover the spatiotemporal patterns of air temperature, LST, and NDVI on the TP across different timescales. These results provide significant insights into the annual and seasonal patterns of vegetation responses and feedback to climate change on the TP. Furthermore, our analysis reveals distinct seasonal trends between NDVI and LST, which can be leveraged to enhance the accuracy of numerical simulations that aim to replicate the relationships between vegetation and climate over the TP.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/rs15092475/s1, Section S1: Detailed procedures of SQMK test; Section S2: Detailed procedures of Sen's slope; Section S3: Overall trends of NDVI, LST, and albedo over the Tibetan Plateau (Figure S1: MODIS-based NDVI, LST, and albedo trends over the Tibetan Plateau at annual and seasonal scales during 2000–2021) [59,60].

**Author Contributions:** Conceptualization, F.W. and Y.M.; methodology, F.W.; manuscript, F.W. and Y.M.; funding acquisition, Y.M. Writing—review and editing, F.W., Y.M., R.D. and C.H.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (2019QZKK0103) and the National Natural Science Foundation of China (42230610 and 91837208).

**Data Availability Statement:** All input data for this research are publicly available: global inventory modeling and mapping studies (GIMMS) NDVI data (https://climatedataguide.ucar.edu/climate-data/ ndvi-normalized-difference-vegetation-index-3rd-generation-nasagfsc-gimms), moderate-resolution imaging spectroradiometer (MODIS) NDVI data (https://lpdaac.usgs.gov/products/mod13a2v061), moderate-resolution imaging spectroradiometer (MODIS) LST data (https://lpdaac.usgs.gov/products/ mod11a2v061), moderate-resolution imaging spectroradiometer (MODIS) albedo data (https://lpdaac. usgs.gov/products/mcd43a3v061), moderate-resolution imaging spectroradiometer (MODIS) Land Cover Type data (https://lpdaac.usgs.gov/products/mcd12q1v061), and surface-near air temperature (https://www.ncdc.noaa.gov/cdo-web).

**Acknowledgments:** The author would like to acknowledge all the members of the Land–Atmosphere Interaction and Its Climatic Effects Group. We also thank reviewers for their valuable comments and the China Scholarship Council (CSC).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
