*2.2. Data Processing*

In this study, alpine meadow and steppe in this region were selected as the research object by using the Vegetation Atlas of China (1:1,000,000). Two NDVI datasets were used, including SPOT NDVI from the Resources and Environment Science and Data Center of the Chinese Academy of Sciences, and MODIS NDVI from the NASA website (MOD13A3 datasets). The spatial resolution of the two types of NDVI is 1 km, with the same time interval of 2000–2018. The maximum NDVI during the growing season (May to September) was used in this study. Temperature and precipitation datasets during the early plant-growing season of 2000–2018 were downloaded from the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn, accessed on 20 July 2019). Considering that the maximum NDVI generally occurs at the peak of the plant-growing season (early August), the climatic conditions in the early growing season (May to July) are critical for vegetation growth [24,25]. Therefore, the precipitation (GSP) and average temperature (GST) of the early growing season were selected, and the spatial resolution was 0.1◦ × 0.1◦. The number of livestock (LN) during 2000–2018 in each county of the study area can be found in the National Tibetan Plateau Data Center and Tibet Statistical Yearbook (www.shujuku.org, accessed on 12 February 2020). The absolute numbers of different animals were converted to standard sheep units (SU) [26]. The vector data, such as administrative divisions used in this study, are from the National Catalog Service for Geographic Information. The population data of the study area were obtained from the National Tibetan Plateau Data Center.

### *2.3. Analysis of NDVI Interannual Change Trend*

A linear regression method was used to detect the interannual changes in the vegetation index [27]. The least square method was used to calculate the regression coefficient, namely, the interannual variation trend of NDVI, as follows:

$$\theta\_{slops} = \frac{\mathbf{n} \times \sum\_{i=1}^{n} \mathbf{i} \times \text{NDVI}\_i - \sum\_{i=1}^{n} \mathbf{i} \sum\_{i=1}^{n} \text{NDVI}\_i}{\mathbf{n} \times \sum\_{i=1}^{n} \mathbf{i}^2 - \left(\sum\_{i=1}^{n} \mathbf{i}\right)^2} \tag{1}$$

where θ*slope* refers to the interannual trend in NDVI, n is the number of years simulated, and NDVI*<sup>i</sup>* is the NDVI value of the *i*th year. We determined the significance of variation (significant change or no significant change) via T-tests to calculate confidence levels.

### *2.4. Correlation Analysis between NDVI and Main Influencing Factors*

Correlation coefficient calculations and significance tests were used to detect the relationship between grassland growth change and climate change in the study area. The partial correlation method was further used to analyze the relative impact of grazing and climate change on NDVI changes. The mathematical expression of the Pearson correlation coefficient is:

$$\mathbf{R\_{X \text{NDVI}}} = \frac{\mathbf{n} \times \sum\_{i=1}^{n} \mathbf{X\_i} \times \text{NDVI}\_i - \sum\_{i=1}^{n} \mathbf{X\_i} \times \sum\_{i=1}^{n} \text{NDVI}\_i}{\sqrt{\mathbf{n} \times \sum\_{i=1}^{n} \mathbf{X\_i}^2 - \left(\sum\_{i=1}^{n} \mathbf{X\_i}\right)^2} \times \sqrt{\mathbf{n} \times \sum\_{i=1}^{n} \text{NDVI}\_i^2 - \left(\sum\_{i=1}^{n} \text{NDVI}\_i\right)^2}} \tag{2}$$

where *R* refers to the correlation coefficient, *n* is the number of years simulated, and X and NDVI are independent and dependent variables, respectively. Based on the NDVI during 2000–2018 and the corresponding temperature and precipitation, this study analyzed the correlation between NDVI and the main climatic factors in the study area, and conducted a significance test. When there are only three variables, a simple correlation coefficient can be used to indirectly calculate the partial correlation coefficient. Setting the three variables as LN, GST, and NDVI, and when LN remains constant, the partial correlation coefficients of GST and NDVI are:

$$\text{R}\_{\text{NDVI GST}} = \frac{\text{R}\_{\text{NDVI GST}} - \text{R}\_{\text{GST LN}} \times \text{R}\_{\text{NDVI LN}}}{\sqrt{(1 - \text{R}\_{\text{GST LN}}^2) \times \left(1 - \text{R}\_{\text{NDVI LN}}^2\right)}}\tag{3}$$

where RNDVI GST refers to the partial correlation coefficient of NDVI and GST when LN is fixed, and where RNDVI SPEI, RNDVI LN, and RGST LN are the correlation coefficients of NDVI and GST, NDVI and LN, and GST and LN, respectively.

### *2.5. RESTREND Analysis*

The residual trend (RESTREND) method is based on the hypothesis that vegetation growth is determined by climate change. After removing the climate influence, the humaninduced vegetation change could thus be identified [28,29]. The trends in NDVI residuals (NDVIH = NDVIP − NDVIA) between the observed (NDVIA) and the predicted values (NDVIP) were calculated at a 0.1◦ × 0.1◦ pixel level. NDVIA is applied to regress against GST and GSP to generate a statistical model which is then used to compute NDVIP at each pixel [30–32]. More details can be found in Li et al. (2012) [30] and Cai et al. (2015) [32]. Table 1 shows the relative effects of climate change (*Slope*NDVIP) and human activities (*Slope*NDVIH) on grassland change (*Slope*NDVIA) in six different scenarios [33]. For example, in the first scenario, an increasing trend in NDVIA and a declining trend in NDVIP result in a decreasing NDVIH, which means that human activity mainly determines the NDVIA increase.


**Table 1.** Evaluation methods of climate change and human activities on grassland dynamics under six possible scenarios.

### *2.6. Statistical Analysis*

The annual variations of NDVI and corresponding climatic factors and the relationships between NDVI and climatic factors were fixed by simple linear correlations. Partial correlations between NDVI and climatic factors of GSP and GST were analyzed at a pixel scale. Since LN data were only available for each county, we applied partial correlation between NDVI and GSP, GST, and LN at a county scale. All statistical analyses were performed using SPSS 19.0 (SPSS Inc., Chicago, IL, USA), and spatial analysis was performed using ArcGIS 10.7 for Desktop. Further, the relative effects of climate change and human activities on NDVI changes were detected by the RESTREND method described above.

### **3. Results**

### *3.1. Interannual Variations in Grassland NDVI during 2000–2018*

Both MODIS NDVI and SPOT NDVI data indicate that the grassland NDVI in the study area increased significantly during 2000–2018 (Figure 2a, MODIS, R2 = 0.28, p = 0.02; SPOT, R2 = 0.33, P = 0.01). The mean NDVI values were approximately 0.23 during the last 19 years, with their maximum in 2013 and minimum in 2002 or 2009. Figure 2b shows that the MODIS NDVI is highly correlated with the SPOT NDVI in the study area (R<sup>2</sup> = 0.99, P = 0.00), indicating that the two datasets can be verified by each other. In view of their high correlation, we only present the results based on MODIS NDVI analysis in the following section.

**Figure 2.** (**a**) Interannual variations in MODIS NDVI and SPOT NDVI and (**b**) their relationships.

### *3.2. Variations in Climatic Factors and the Number of Livestock and Their Relationships with NDVI*

In this area, GST (Figure 3a, R2 = 0.02, P = 0.55) and GSP (Figure 3b, R<sup>2</sup> = 0.03, P = 0.52) did not vary during 2000–2018. However, LN reflected a significant decreasing trend (Figure 3c, R2 = 0.56, P = 0.00), especially after 2008.

**Figure 3.** Interannual variations in (**a**) GST, (**b**) GSP, and (**c**) LN.

Pearson correlation analysis showed that NDVI was not correlated with GST (Figure 4a, R2 = 0.02, P = 0.57) but was positively related to GSP (Figure 4b, R<sup>2</sup> = 0.35, P = 0.01) and negatively correlated with LN (Figure 4c, R2 = 0.23, P = 0.04).

**Figure 4.** The relationships between NDVI and (**a**) GST, (**b**) GSP, and (**c**) LN.

### *3.3. Relative Effects of Climate Factors and Grazing on Grassland NDVI*

At the pixel level, NDVI was more closely related with GSP than with GST. There were 22.31% pixels showing significantly positive partial correlations between NDVI and GSP (Figure 5a), but only 5.42% reflecting significantly positive correlations between NDVI and GST (Figure 5b). Additionally, the partial correlation coefficients between NDVI and GSP showed substantial spatial variations, and those area affected by GSP were mainly distributed in the center part, i.e., Saga and Ngamring.

At the regional level, the partial correlation analysis showed that NDVI variations for most counties in this area were mainly affected by the climatic factor of precipitation, and only that in Tingri County were mediated by both climate and grazing (Figure 6). Among the six counties affected by climate factors, the NDVI for Dinggye County was mainly controlled by GST, while the other five were generally affected by GSP.

**Figure 5.** Spatial distribution of the partial correlation coefficients between NDVI and (**a**) GSP, and (**b**) GST. Different colors indicate the classes of partial correlation coefficients. Different symbols denote significant levels: \* *p* < 0.05; # 0.05 < *p* < 0.10; ns, not significant. The map was created using ESRI ArcMap 10.7 software; the topographic base of the map was created with National Catalogue Service For Geographic Information data (https://www.webmap.cn, figure number: GS(2016)2556).

**Figure 6.** Main drivers of NDVI change in the seven counties of the study area. Different symbols indicate the different significant levels (\* *p* < 0.05; # 0.05 < *p* < 0.10) according to the partial correlation coefficients. The map was created using ESRI ArcMap 10.7 software; the topographic base of the map was created with National Catalogue Service For Geographic Information data (https: //www.webmap.cn, figure number: GS(2016)2556).

The RESTREND analysis, as shown by Figure 7, further indicated that the area showing increasing NDVI accounted for 63.22% of the study area, in which 32.17% was due to climate change, 24.80% was due to human activities, and 6.25% was due to both climate change and human activities. The reason for the NDVI increase by human activities might be mainly due to the implementation of an ecological restoration project (Grazing Withdrawal Program). For the area (36.78%) depicting decreasing NDVI, 21.62% was due to human activities, 12.30% due to climate change, and 2.87% due to both climate change and human activities. The reason for the NDVI decrease by human activities could be mainly due to overgrazing (e.g., Tingri County).

**Figure 7.** Relative effects of climate change and human activities on NDVI changes in the Mt. Qomolangma Nature Reserve and adjacent regions during 2000–2018. HI denotes human-activityinduced NDVI increase; HPI denotes both climate-change and human-activity-induced NDVI increase; PI denotes climate-change-induced NDVI increase; HD denotes NDVI decrease induced by human activities; HPD denotes NDVI decrease induced by both climate change and human activities; and PD denotes NDVI decrease induced by climate change. The map was created using ESRI ArcMap 10.7 software; the topographic base of the map was created with National Catalogue Service For Geographic Information data (https://www.webmap.cn, figure number: GS(2016)2556).

### **4. Discussion**

The Tibetan Plateau has an average altitude of more than 4000 m and is characterized as a typical alpine ecosystem. Therefore, low temperature is generally considered to be a key factor restricting plant growth in this area [6,10]. However, in arid western Tibet, precipitation is considered to be a more important limiting factor than temperature [34]. Synthesizing the results of Pearson correlation and partial correlation analysis, we found that the variations in grassland NDVI in the Mt. Qomolangma National Nature Reserve and adjacent regions are mainly affected by climatic factors, especially precipitation. This is consistent with existing studies. For example, based on an analysis of the net primary productivity (NPP) of grassland vegetation on the Tibetan Plateau during 2000–2015, He et al. (2019) found that precipitation was the first limiting factor of NPP [35]. This is associated with the effects of rain shadows in the Himalayas. At the same time, the impact of climate change on grassland vegetation is largely related to grazing intensity. In areas with heavy grazing intensity (larger LN), grassland NDVI is more significantly affected by grazing, and when grazing intensity exceeds a certain threshold, the negative effect of grazing on grassland is likely to offset or even exceed the positive effect of climate change on the growth of alpine grassland [36]. In general, human activities have a relatively weak impact on the grassland ecosystems in the study area; therefore, the variation in NDVI is mainly controlled by climate.

However, unlike the other six counties, the variation in grassland NDVI in Tingri County was controlled by both climate (precipitation) and grazing, which indicated that the grazing pressure was heavy in Tingri County to some extent. Figure 8a further illustrated that the population density of Tingri County was the highest among the seven counties, and the grazing intensity was the second highest. Synthesizing the population and grazing intensity data of the seven counties in the study area during the past 20 years, we found that both data were significantly correlated with each other (Figure 8b), i.e., the variation

in population density at the county level largely reflected that of the grazing intensity. In addition, the Mt. Qomolangma scenic spot is located in Tingri County. According to the data provided by the Mt. Qomolangma National Nature Reserve Administration, an average of 82,300 tourists was received each year during the past 10 years (2011–2019). This number increased to approximately 119,000 during the past 3 years (2017–2019), which is approximately twice the resident population of Tingri County. A large number of tourists would undoubtedly stimulate the demands for livestock products, which might indirectly cause an increase in grazing intensity and aggravate the negative impact on grassland.

**Figure 8.** (**a**) Comparisons of the grazing intensity and population density between 7 counties and (**b**) their relationships during 2000–2018. Different lower-case and upper-case letters indicate significant differences in grazing intensity and population density between counties, respectively, by using the Tukey's HSD tests at α = 0.05.

It is worth noting that the grazing intensity of Dinggye County ranks first among the seven counties in the study area (Figure 8a). Theoretically, the grazing intensity should be closely related to the changes in NDVI, but partial correlation analysis showed that the grazing effect on the grassland NDVI of the county was not significant. This may be related to the vegetation degradation caused by overgrazing in the county; that is, under heavy grazing conditions, the relationship between NDVI and grazing intensity may be nonlinear [37,38]. In addition, it can be inferred from Figure 9 that NDVI showed a significant negative correlation with grazing intensity after excluding the data in 2008 and 2011, when the precipitation was abundant or the temperature was very low, indicating that the strong fluctuation of the hydrothermal combination factor in different years may play an important role in regulating the relationship between grazing intensity and NDVI. The significantly negative partial correlation between GST and NDVI in Dinggye County further indicates that the future warming climate will not be advantageous for plant growth, which is assumed to associate with increased evapotranspiration and droughts [39,40].

In this study, we found that the NDVI showed an overall increasing trend during 2000– 2018, and GSP was the main limiting factor. However, the GSP did not show a significant increasing pattern in the past 19 years, indicating that its interannual fluctuations rather than the changing trend determined the correlation between NDVI and precipitation. This result is consistent with Zhang et al. (2013), who found that the decline in NDVI in the Kosi River Basin in the central Himalayas (involving Tingri, Dinggye, and Nyalam counties) during 1982–1994 and 1994–2000 was related to a sudden drop in precipitation [3]. Therefore, concerning the reasons for vegetation changes, we should pay more attention to the impact of extreme climate events such as extreme precipitation years [41,42].

**Figure 9.** Relationship between NDVI and grazing intensity in Dinggye County.
