**1. 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

**Citation:** Wang, F.; Ma, Y.; Darvishzadeh, R.; Han, C. Annual and Seasonal Trends of Vegetation Responses and Feedback to Temperature on the Tibetan Plateau since the 1980s. *Remote Sens.* **2023**, *15*, 2475. https://doi.org/10.3390/ rs15092475

Academic Editor: Conghe Song

Received: 6 April 2023 Revised: 2 May 2023 Accepted: 5 May 2023 Published: 8 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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

### **2. Study Area and Datasets**

*2.1. Study Area*

With an area of 257 × <sup>10</sup><sup>4</sup> km<sup>2</sup> 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.

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

### *2.2. 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.



### 2.2.1. Meteorological Station Data

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.

### 2.2.2. 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–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].

### *2.3. 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–41], the Sequential Mann–Kendall (SQMK) [32–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.

### **3. Results**

*3.1. Air Temperature and Vegetation Trends during 1982–2013 at Annual and Seasonal Scales* 3.1.1. Annual Trends in Air Temperature and Vegetation

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.

**Figure 2.** *Cont*.

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

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

**Figure 3.** *Cont*.

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

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), particularly in the western TP.

**Figure 4.** *Cont*.

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

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,f), which is similar to the vegetation change pattern in spring.

**Figure 5.** *Cont*.

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

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

**Figure 6.** *Cont*.

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

### 3.1.3. 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.
