Annual and Seasonal Trends of Vegetation Responses and Feedback to Temperature on the Tibetan Plateau since the 1980s

: The vegetation–temperature relationship is crucial in understanding land–atmosphere interactions on the Tibetan Plateau. Although many studies have investigated the connections between vegetation and climate variables in this region using remote sensing technology, there remain notable gaps in our understanding of vegetation–temperature interactions over different timescales. Here, we combined site-level air temperature observations, information from the global inventory modeling and mapping studies (GIMMS) dataset, and moderate-resolution imaging spectroradiometer (MODIS) products to analyze the spatial and temporal patterns of air temperature, vegetation, and land surface temperature (LST) on the Tibetan Plateau at annual and seasonal scales. We achieved these spatiotemporal patterns by using Sen’s slope, sequential Mann–Kendall tests, and partial correlation analysis. The timescale differences of vegetation-induced LST were subsequently discussed. Our results demonstrate that a breakpoint of air temperature change occurred on the Tibetan Plateau during 1994–1998, dividing the study period (1982–2013) into two phases. A more signiﬁcant greening response of NDVI occurred in the spring and autumn with earlier breakpoints and a more sensitive NDVI response occurred in recent warming phase. Both MODIS and GIMMS data showed a common increase in the normalized difference vegetation index (NDVI) on the Tibetan Plateau for all timescales, while the former had a larger greening area since 2000. The most prominent trends in NDVI and LST were identiﬁed in spring and autumn, respectively, and the largest areas of signiﬁcant variation in NDVI and LST mostly occurred in winter and autumn, respectively. The partial correlation analysis revealed a signiﬁcant negative impact of NDVI on LST during the annual scale and autumn, and it had a signiﬁcant positive impact during spring. Our ﬁndings improve the general understanding of vegetation–climate relationships at annual and seasonal scales.


Introduction
Vegetation is an important part of the global terrestrial ecosystem and can significantly impact global physical energy cycles, carbon balance regulations, and regional climate [1,2].Global vegetation cover generally increases with the warming of the Earth's climate [3], and the surface albedo changes caused by vegetation change affect both the net amount of solar radiation absorbed by the earth surface and the evapotranspiration rate.Consequently, vegetation-related evapotranspiration and albedo affect vegetation-temperature relationships [4,5].The Tibetan Plateau (TP) is the highest plateau on earth and possesses a complex topography, which generates a unique relationship between vegetation and temperature that is extremely sensitive to global climate change [6].Investigating the responses and feedbacks of vegetation to the TP's temperature is crucial for understanding land-atmosphere interactions.
Numerous researchers have studied the vegetation-temperature relationships on the TP through different approaches, including field observations [7,8], numerical simulations, and remote sensing observations [9,10].However, given the challenging topography and extreme climate of the TP [11], fieldwork and in situ analysis can be logistically difficult to perform.It is not easy to consistently assess vegetation-temperature relationships at different timescales.As such, many studies have employed numerical models to simulate and analyze their relationships [12], and significant research progress has been made using these techniques.However, since models often have considerable uncertainties and coarse spatial resolution, it is unfeasible to numerically simulate the vegetation-temperature relationships at different timescales for long periods [13].Satellite Earth observation is an effective instrument for monitoring bio-geophysical variables of vegetation at large scales [14].The vegetation index and albedo derived from satellite observations can capture vegetation greening and browning and surface albedo change.These changes alter surface bio-geophysical properties and near-surface aerodynamics, leading to an effect on local temperature through bio-geophysical feedbacks [15][16][17].As satellite remote sensing technology has developed, the products of widespread vegetation-temperature interactions that accumulate on the TP can be measured remotely, allowing a quantitative investigation of vegetation coverage at large scales and over extended periods [6,18].Consequently, remote sensing observations are becoming a valuable instrument for investigating vegetation and climate interactions.
Using the available remote sensing data, previous papers report the trends and the spatial variability of vegetation, surface albedo, and land surface temperature (LST) on the TP at annual scales.The satellite-derived normalized difference vegetation index (NDVI) indicates the status of plant growth.This indicator can be quantified by the difference between the near-infrared (representing vegetation reflection) and red bands (representing vegetation absorption) [19].Global inventory modeling and mapping studies (GIMMS) data suggest an increasing NDVI trend on the TP during 1985-1999, which is likely due to the shift from an arid to a warm-humid climate and the reduction in human activities in this region [20].Similar studies have found an increasing NDVI trend during the latter two decades of the 20th century [21], which has been confirmed by researchers since the year 2000.Statistical analysis of SPOT (Satellite Pour l' Observation de la Terre) data collected during 1999-2014 shows an overall increasing trend in NDVI on the TP coupled with a moderate increase in air temperature [6].Studies based on MODIS data collected over the TP during 2001-2019 also support this, as they show an increasing trend in NDVI and LST but a decreasing trend in albedo during that period.Increased forest coverage and decreased snow coverage are considered to be the dominant factors that drove these changes [22].Further studies have shown that increased snowfall induced an increase in albedo on the southwestern TP due to anomalous water vapor transport [23].Indeed, the LST warmed at a significantly faster rate than the air temperature, with the annual temperature increase during 1987-2008 in the former showing 0.78 ± 0.0631 K/decade, but in the latter, it showed 0.275 ± 0.0216 K/decade [24].LST on the TP is influenced by various factors such as elevation, surface radiation, subsurface temperature, and surface properties [25][26][27].Nevertheless, vegetation and albedo are becoming the hot spots for LST warming studies.This can be attributed to the significant linear relationship between vegetation and LST [28], the direct contribution of surface albedo to LST [25], and satellite advantages [29].Despite different spatial patterns being inferred from different datasets, they show consistent support for the overall observed trends and patterns of vegetation and temperature on the TP.
Most of the studies mentioned above focused on the annual scale, which leaves two important issues unresolved: (1) no systematic analysis has been performed to analyze the differences in both the warming phases (well known as a period of increasing temperature) and the vegetation response to temperature across different timescales over the TP; and (2) the spatial-temporal pattern of vegetation and its effects on LST across the TP remains unclear.We addressed these problems by combining site-level observations with GIMMS data to (1) analyze spatial and temporal patterns at the annual and seasonal scales of air temperature and NDVI during different warming phases in 1982-2013 and (2) examine the annual and seasonal scales of NDVI and LST changes during 2000-2021 using MODIS data products and thus discuss the possible impacts of vegetation on LST at the annual and seasonal scales.Our results have revealed systematic annual and seasonal characteristics of air temperature, LST, and vegetation changes on the TP, and this can be used to form a better understanding of land-atmosphere interaction patterns across the TP.

Study Area
With an area of 257 × 10 4 km 2 and an average altitude above 4000 m, the TP is the largest and highest plateau in the world, and it exhibits a complex topography and diverse underlying surface conditions (Figure 1).The climate and ecological environment on the TP have both changed dramatically since the middle of the 20th century due to intensified human activity [30].Significant warm and wet trends on the TP have occurred since the 1960s as documented by temperature and precipitation data [31].Due to its distinctive geographic location and susceptibility to climate change, the TP serves as a natural laboratory for investigating the intricate relationship between vegetation and temperature.
important issues unresolved: (1) no systematic analysis has been performed to analyze the differences in both the warming phases (well known as a period of increasing temperature) and the vegetation response to temperature across different timescales over the TP; and (2) the spatial-temporal pattern of vegetation and its effects on LST across the TP remains unclear.We addressed these problems by combining site-level observations with GIMMS data to (1) analyze spatial and temporal patterns at the annual and seasonal scales of air temperature and NDVI during different warming phases in 1982-2013 and (2) examine the annual and seasonal scales of NDVI and LST changes during 2000-2021 using MODIS data products and thus discuss the possible impacts of vegetation on LST at the annual and seasonal scales.Our results have revealed systematic annual and seasonal characteristics of air temperature, LST, and vegetation changes on the TP, and this can be used to form a better understanding of land-atmosphere interaction patterns across the TP.

Study Area
With an area of 257 × 10 4 km 2 and an average altitude above 4000 m, the TP is the largest and highest plateau in the world, and it exhibits a complex topography and diverse underlying surface conditions (Figure 1).The climate and ecological environment on the TP have both changed dramatically since the middle of the 20th century due to intensified human activity [30].Significant warm and wet trends on the TP have occurred since the 1960s as documented by temperature and precipitation data [31].Due to its distinctive geographic location and susceptibility to climate change, the TP serves as a natural laboratory for investigating the intricate relationship between vegetation and temperature.

Data Sources and Processing
The data used in this study were obtained from meteorological station observations (Figure 1) (near-surface air temperature) and remote sensing (Table 1).Remote sensing data comprised NDVI products derived from GIMMS and MODIS, albedo, land cover type, and LST products derived from MODIS.In this study, winter, spring, summer, and autumn are defined as extending from December to the following February (DJF), March to May (MAM), June to August (JJA), and September to November (SON), respectively.Air temperature data were collected from the National Climate Data Center, which is affiliated with the National Oceanic and Atmosphere Administration (NOAA) (Figure 1).These data were used to produce the average air temperature at different timescales from 1982 to 2013.First, we acquired 44 high-quality meteorological stations considering the geographical scope of the TP and the requirement of high-density temperature data.Second, we collected the daily mean temperature data for each station during 1982-2013 and replaced the missing values for a few of stations using the average air temperature of the neighboring days [32,33].Finally, the mean temperature of annual and seasonal scales was calculated and used to detect trends and breakpoints in different warming phases of the TP.

Remote Sensing Products
The remote sensing data used in this study included the GIMMS and MODIS datasets.The GIMMS NDVI product (1982-2013) was derived from the AVHRR sensor of the National Oceanic and Atmospheric Administration (http://www.ncdc.noaa.gov/cdo-web,accessed on 3 October 2022), and it has a temporal resolution of 15 days and a spatial resolution of 8 km.The initial version of this dataset is not ideal for capturing vegetation dynamics, and therefore, the latest version mitigates this problem by correcting sensors, aerosols, and view geometry [34][35][36].A maximum value composite procedure was used to remove some sources of interference, such as clouds, the atmosphere, and variation in solar altitude angle; after that, annual and seasonal NDVI values were obtained [37].
Compared to the AVHRR instrument, the updated MODIS instrument has a better sensitivity to chlorophyll with higher spatial resolution.We employed MODIS datasets of 2000-2021 including NDVI (MOD13A2), LST (MOD11A1), albedo (MCD43A3), and land cover types (MCD12Q1).Specifically, MOD13A2 and MOD11A1 have a 1 km spatial resolution and a 16-day temporal resolution.MOD11A1 and MCD43A3 were used in this study to better understand the relationships between LST and vegetation change.MCD12Q1 (version 6.1), providing annual land cover types (2001-2021), was also obtained to analyze the effect of land cover change on NDVI change trends.Additionally, all of the MODIS datasets were aggregated to 1 km to ensure the consistent spatial resolution of these datasets.We adopted different strategies in processing the MODISderived surface parameters by referring to previous studies.For each timescale, we processed the NDVI data using the maximum value composite method [37], and we processed the LST and albedo data using the mean value composite method [29,38].

Methods
Our analyses were performed in three major steps.Firstly, to analyze the NDVI responses at different warming phases, we performed breakpoint detection of air temperature data.Considering previous artificial temporal segmentation on NDVI variation [17,[39][40][41], the Sequential Mann-Kendall (SQMK) [32][33][34]42] was employed to estimate the breakpoints of temperature change, and it was used as a proxy for classifying the different phases of NDVI responses.Please see the Supporting Information S1 for the detailed introduction to the SQMK test.
Secondly, we analyzed trends and performed a significance analysis [43,44] on key surface parameters.Sen's slope estimator, a nonparametric test [45], was employed to determine trends in NDVI, LST, and albedo for different warming phases of all timescales.Please see the Supporting Information S2 for the detailed introduction to the Sen's slope.
Thirdly, to identify the impacts of the NDVI on LST, both the detrending method and partial correlation analysis were performed among the NDVI, albedo, and LST.The purpose of using the detrending method was to eliminate any spurious correlations among the three parameters that may have been caused by temporal variations.As for the partial correlation analysis, it was used to better quantify the individual impact of NDVI or albedo on LST.The specific procedures included two aspects: (1) Using the first-difference detrending method (i.e., the difference of values in one year to the previous year), we examined and filtered the temporal trends of NDVI and albedo.(2) For partial correlation analysis, the partial correlation coefficient is an assessment of the net correlation between a single factor and the target value, provided that the impact of other factors is fixed or deducted.Considering the significant co-impact of vegetation and albedo on LST, the partial correlation coefficient is a good indicator for analyzing the relationship between them.The air temperature and vegetation coverage on the TP generally increased on an annual scale (Figure 2a), while vegetation changes during each warming phase showed significant spatial and temporal variation.A significant abrupt change (p < 0.05) occurred in 1996 for air temperature trends on the TP (Figure 2b), and the warming trend during 1996-2013 (0.043 • C/year) was notably higher than that during 1982-1996 (0.042 • C/year), suggesting that the warming rate on the TP accelerated after 1996.The significant result (p < 0.05) revealed that the annual NDVI is generally increasing (Figure 2a) and the clustered NDVI increase and decrease in the second warming phase are greater than that in the first warming phase (Figure 2d,f).Specifically, the NDVI greened more than it browned during the first warming phase (Figure 2d), and there was clustered NDVI greening and browning in the eastern and western plateau (Figure 2f), respectively.

Air Temperature and Vegetation
Remote Sens. 2023, 15, x FOR PEER REVIEW 6 of

Seasonal Trends in Air Temperature and Vegetation
The air temperature and vegetation coverage on the TP during the spring seaso generally showed an increasing trend throughout the study period (Figure 3a).The brea point and warming rates were very different from those of the annual scale.The brea point for spring's air temperature trend occurred in 1994, with a slope of 0.040 °C/year the first warming phase and 0.034 °C/year in the second warming phase (Figure 3b), ind cating that warming on the TP during the spring has slowed down.The NDVI showed general increase over time, although the trend was more significant during the secon warming phase (Figure 3c,e).Furthermore, the NDVI showed a significant increasin trend (p < 0.05) in the eastern TP and a significant decreasing trend in the northweste part (Figure 3e).The significance statistics (p < 0.05) suggested that the clustered NDV changes occurred in the second phase rather than in the first phase (Figure 3d,f).For e ample, TP NDVI decreased in the western part and increased in the eastern part (Figu 3f).

Seasonal Trends in Air Temperature and Vegetation
The air temperature and vegetation coverage on the TP during the spring seasons generally showed an increasing trend throughout the study period (Figure 3a).The breakpoint and warming rates were very different from those of the annual scale.The breakpoint for spring's air temperature trend occurred in 1994, with a slope of 0.040 • C/year in the first warming phase and 0.034 • C/year in the second warming phase (Figure 3b), indicating that warming on the TP during the spring has slowed down.The NDVI showed a general increase over time, although the trend was more significant during the second warming phase (Figure 3c,e).Furthermore, the NDVI showed a significant increasing trend (p < 0.05) in the eastern TP and a significant decreasing trend in the northwestern part (Figure 3e).The significance statistics (p < 0.05) suggested that the clustered NDVI changes occurred in the second phase rather than in the first phase (Figure 3d,f).For example, TP NDVI decreased in the western part and increased in the eastern part (Figure 3f).For summer, the TP exhibited both the largest NDVI value and the most significan warming trend (Figure 4a, p < 0.05), as did the largest pre-and post-phase difference i air temperature (Figure 4b).The breakpoint for the summer temperature on the TP oc curred in 1998, which is slightly later than for the spring and annual scales.The trend fo the second warming phase reached a rate of 0.047 °C/year, which is notably higher tha that of the first warming phase.This implied that the warming rate in summer on the T is much greater than that in other timescales.The significance analysis (p < 0.05) indicate that the NDVI increase is greater in the first warming phase than that in the second warm ing phase (Figure 4d,f), particularly in the western TP.For summer, the TP exhibited both the largest NDVI value and the most significant warming trend (Figure 4a, p < 0.05), as did the largest pre-and post-phase difference in air temperature (Figure 4b).The breakpoint for the summer temperature on the TP occurred in 1998, which is slightly later than for the spring and annual scales.The trend for the second warming phase reached a rate of 0.047 • C/year, which is notably higher than that of the first warming phase.This implied that the warming rate in summer on the TP is much greater than that in other timescales.The significance analysis (p < 0.05) indicated that the NDVI increase is greater in the first warming phase than that in the second warming phase (Figure 4d,f  Autumn vegetation trends and warming both showed a lower value compared to th summer (Figure 5a).The breakpoint of temperature during the autumn occurred in 199 (Figure 5b), similar to the spring.The first and second phases during autumn recorde trends of 0.029 °C/year and 0.018 °C/year, respectively.As such, the trend for the secon phase was significantly weaker than that of the first phase, demonstrating that the warm ing rate on the TP becomes slow during autumn.The NDVI showed a non-significan fluctuation in general, and the significance statistics (p < 0.05) showed that the increase o decrease in the TP's NDVI is less significant in the first phase than that in the second phas (Figure 5d,f), which is similar to the vegetation change pattern in spring.Autumn vegetation trends and warming both showed a lower value compared to the summer (Figure 5a).The breakpoint of temperature during the autumn occurred in 1994 (Figure 5b), similar to the spring.The first and second phases during autumn recorded trends of 0.029 • C/year and 0.018 • C/year, respectively.As such, the trend for the second phase was significantly weaker than that of the first phase, demonstrating that the warming rate on the TP becomes slow during autumn.The NDVI showed a non-significant fluctuation in general, and the significance statistics (p < 0.05) showed that the increase or decrease in the TP's NDVI is less significant in the first phase than that in the second phase (Figure 5d  The vegetation and air temperature changes on the TP documented during the winte were similar to those during the autumn, although they showed much smaller magn tudes than other timescales (Figure 6a,b).The breakpoint in the winter temperature be tween warming phases on the TP occurred in 1998; the trends of the first and secon warming phases were 0.017 °C/year and 0.021 °C/year, respectively, indicating that th winter warming on the TP is accelerating (Figure 6b).Winter changes in the NDVI value were generally small.A fragmented NDVI increase occurred in the eastern TP of the fir phase, but a massive NDVI decrease occurred in the second phase, which was not foun for any other timescales (Figure 6d,f).The vegetation and air temperature changes on the TP documented during the winter were similar to those during the autumn, although they showed much smaller magnitudes than other timescales (Figure 6a,b).The breakpoint in the winter temperature between warming phases on the TP occurred in 1998; the trends of the first and second warming phases were 0.017 • C/year and 0.021 • C/year, respectively, indicating that the winter warming on the TP is accelerating (Figure 6b).Winter changes in the NDVI values were generally small.A fragmented NDVI increase occurred in the eastern TP of the first phase, but a massive NDVI decrease occurred in the second phase, which was not found for any other timescales (Figure 6d,f  Our study identified consistent breakpoints in air temperature for spring and au tumn as well as for summer and winter.We observed that the first warming phase showe a smaller trend than the second warming phase on most timescales.This suggests th while warming continues on the TP, the later warming trend is slightly decreasing.
We also investigated the spatial trends of NDVI at different warming phases.Durin the first warming phase, we found that the area of NDVI greening is smaller in spring an autumn with earlier breakpoints compared to summer and winter with later breakpoint However, in the second warming phase, we observed that the area of NDVI increase larger in spring and autumn with earlier breakpoints compared to summer and winte with later breakpoints.This indicates that the timescale with the earlier breakpoints e hibits a greater NDVI increase in the second warming phase.For instance, the air tempe ature breakpoint at the annual scale earlier than summer and winter showed a signif cantly larger area of NDVI increase in the second warming phase.
Furthermore, we found that in the relatively weaker second warming phase, there more NDVI decrease at all timescales across the sparsely vegetated northwestern TP.Th may be attributed to two factors: (1) limitations in the ability of the AVHRR sensor t capture detailed NDVI changes in sparsely vegetated areas, and (2) an increased veget tion sensitivity to air temperature during the second warming phase.

The Spatial and Temporal Responses of NDVI to Air Temperature
Our study identified consistent breakpoints in air temperature for spring and autumn as well as for summer and winter.We observed that the first warming phase showed a smaller trend than the second warming phase on most timescales.This suggests that while warming continues on the TP, the later warming trend is slightly decreasing.
We also investigated the spatial trends of NDVI at different warming phases.During the first warming phase, we found that the area of NDVI greening is smaller in spring and autumn with earlier breakpoints compared to summer and winter with later breakpoints.However, in the second warming phase, we observed that the area of NDVI increase is larger in spring and autumn with earlier breakpoints compared to summer and winter with later breakpoints.This indicates that the timescale with the earlier breakpoints exhibits a greater NDVI increase in the second warming phase.For instance, the air temperature breakpoint at the annual scale earlier than summer and winter showed a significantly larger area of NDVI increase in the second warming phase.
Furthermore, we found that in the relatively weaker second warming phase, there is more NDVI decrease at all timescales across the sparsely vegetated northwestern TP.This may be attributed to two factors: (1) limitations in the ability of the AVHRR sensor to capture detailed NDVI changes in sparsely vegetated areas, and (2) an increased vegetation sensitivity to air temperature during the second warming phase.

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).8).These NDVI data showed that the eastern TP records the largest

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).

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.

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].
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).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.

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.

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

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.

Figure 1 .
Figure 1.Location of meteorological stations on the Tibetan Plateau.

Figure 1 .
Figure 1.Location of meteorological stations on the Tibetan Plateau.
Trends during 1982-2013 at Annual and Seasonal Scales 3.1.1.Annual Trends in Air Temperature and Vegetation

Figure 2 .
Figure 2. Vegetation and air temperature trends on the annual scale.GIMMS-based NDVI and a temperature trends during 1982-2013 (a); air temperature breakpoint detection based on the M method (b).In (b), the gray dashed line denotes the detected temperature breakpoint, and UF an UB refer to the statistics of forward and backward sequence, respectively; GIMMS-based ND trends (c) and regions of significant change (d) during the first warming phase; GIMMS-based ND trends (e) and regions of significant change (f) during the second warming phase.

Figure 2 .
Figure 2. Vegetation and air temperature trends on the annual scale.GIMMS-based NDVI and air temperature trends during 1982-2013 (a); air temperature breakpoint detection based on the MK method (b).In (b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refer to the statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends (c) and regions of significant change (d) during the first warming phase; GIMMS-based NDVI trends (e) and regions of significant change (f) during the second warming phase.

Figure 3 .
Figure 3. Vegetation and air temperature trends in spring.GIMMS-based NDVI and air temperatur trends during 1982-2013 (a); air temperature breakpoint detection based on the MK method (b).I (b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refer to th statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends (c) and re gions of significant change (d) during the first warming phase.GIMMS-based NDVI trends (e) an regions of significant change (f) during the second warming phase.

Figure 3 .
Figure 3. Vegetation and air temperature trends in spring.GIMMS-based NDVI and air temperature trends during 1982-2013 (a); air temperature breakpoint detection based on the MK method (b).In (b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refer to the statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends (c) and regions of significant change (d) during the first warming phase.GIMMS-based NDVI trends (e) and regions of significant change (f) during the second warming phase.

Figure 4 .
Figure 4. Vegetation and air temperature trends in summer.GIMMS-based NDVI and air tempera ture trends during 1982-2013 (a); air temperature breakpoint detection based on the MK metho (b).In (b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refe to the statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends (c and regions of significant change (d) during the first warming phase.GIMMS-based NDVI trend (e) and regions of significant change (f) during the second warming phase.

Figure 4 .
Figure 4. Vegetation and air temperature trends in summer.GIMMS-based NDVI and air temperature trends during 1982-2013 (a); air temperature breakpoint detection based on the MK method (b).In (b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refer to the statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends (c) and regions of significant change (d) during the first warming phase.GIMMS-based NDVI trends (e) and regions of significant change (f) during the second warming phase.

Figure 5 .
Figure 5. Vegetation and air temperature trends in autumn.GIMMS-based NDVI and air temper ture trends during 1982-2013 (a); air temperature breakpoint detection based on the MK metho (b).In (b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refe to the statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends ( and regions of significant change (d) during the first warming phase.GIMMS-based NDVI trend (e) and regions of significant change (f) during the second warming phase.

Figure 5 .
Figure 5. Vegetation and air temperature trends in autumn.GIMMS-based NDVI and air temperature trends during 1982-2013 (a); air temperature breakpoint detection based on the MK method (b).In (b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refer to the statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends (c) and regions of significant change (d) during the first warming phase.GIMMS-based NDVI trends (e) and regions of significant change (f) during the second warming phase.

Figure 6 .
Figure 6.Vegetation and air temperature trends in winter.GIMMS-based NDVI and air temperatu trends during 1982-2013 (a); air temperature breakpoint detection based on the MK method (b).(b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refer to th statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends (c) and r gions of significant change (d) during the first warming phase.GIMMS-based NDVI trends (e) an regions of significant change (f) during the second warming phase.

Figure 6 .
Figure 6.Vegetation and air temperature trends in winter.GIMMS-based NDVI and air temperature trends during 1982-2013 (a); air temperature breakpoint detection based on the MK method (b).In (b), the gray dashed line denotes the detected temperature breakpoint, and UF and UB refer to the statistics of forward and backward sequence, respectively; GIMMS-based NDVI trends (c) and regions of significant change (d) during the first warming phase.GIMMS-based NDVI trends (e) and regions of significant change (f) during the second warming phase.

22 Figure 7 .
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 (Figure8).These NDVI data showed that the eastern TP records the largest

Figure 7 .
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.

Table 1 .
Data sources used in this study.

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).