Vegetation Changes and Their Response to Global Change Based on NDVI in the Koshi River Basin of Central Himalayas Since 2000

Abstract

Vegetation forms a main component of the terrestrial biosphere owing to its crucial role in land cover and climate change, which has been of wide concern for experts and scholars. In this study, we used MODIS (moderate-resolution imaging spectroradiometer) NDVI (Normalized Difference Vegetation Index) data, land cover data, meteorological data, and DEM (Digital Elevation Model) data to do vegetation change and its relationship with climate change. First, we investigated the spatio-temporal patterns and variations of vegetation activity in the Koshi River Basin (KRB) in the central Himalayas from 2000 to 2018. Then, we combined NDVI change with climate factors using the linear method to examine their relationship, after that we used the literature review method to explore the influence of human activities to vegetation change. At the regional scale, the NDVI_GS (Growth season NDVI) significantly increased in the KRB in 2000–2018, with significant greening over croplands in KRB in India. Further, the croplands and forest in the KRB in Nepal were mainly influenced by human interference. For example, improvements in agricultural fertilization and irrigation facilities as well as the success of the community forestry program in the KRB in Nepal increased the NDVI_GS of the local forest. Climate also had a certain impact on the increase in NDVI_GS. A significant negative correlation was observed between NDVI_GS trend and the annual minimum temperature trend (TMN) in the KRB in India, but an insignificant positive correlation was noted between it and the total annual precipitation trend (PRE). NDVI_GS significantly decreased over a small area, mainly around Kathmandu, due to urbanization. Increases in NDVI_GS in the KRB have thus been mainly affected by human activities, and climate change has helped increase it to a certain extent.

1. Introduction

Vegetation is the main component of the terrestrial ecosystem, and plays a critical role in the global carbon, water, and energy cycles [1,2]. The NDVI index (Normalized Difference Vegetation Index), a reliable indicator of vegetation change, at both regional and global scales [3], is often used as representative of the status of terrestrial vegetation coverage and growth, because it is closely correlated with the fraction of photosynthetically active radiation absorbed by plant canopies [4], LAI (leaf area index) [5,6], potential photosynthesis [7], and biomass, which is one of the best indicators for quantitatively characterizing vegetation dynamic changes [8]. Changes in vegetation productivity regulate net carbon uptake, and alter the biophysical properties of the land surface and energy transfer, which ultimately leads to local and global changes in climate [9,10]. Therefore, the dynamics of vegetation and its response to climate change have received extensive attention in research.

Studies on this subject have shown that global vegetation has changed in recent decades over various climate zones, vegetation types, and soil types [1]. These changes include the greening in Europe [1], China [11], and India [12], and the browning over northern North America [2] and Southeast Asia [13]. The Himalayan mountain system is among the most dynamic in the world, and is particularly sensitive to climatic change. It features large altitude drops and various types of vegetation. The significant trend in Himalayan vegetation has been one of greening in recent decades, with the distribution exhibiting a strong dependence on elevation and the ecoregion. Greening is most dominant at low and middle elevations, whereas a majority of the browning has occurred at higher elevation, with the eastern high Himalayas browning more prominently than the western high Himalayas [14]. In the eastern part of the Himalayas, prominent browning has occurred in the vegetation at high altitudes. The main greening of dry land and irrigated agricultural areas is determined based on the classification of land cover and land use [14]. Zhang et al. [15] analyzed the spatio-temporal patterns of vegetation growth from 1982 to 2006 in the Koshi River Basin (KRB in China and in Nepal) based on GIMMS (Global Inventory Modeling and Mapping Studies) data at a resolution of 8 km, and used vegetation data from SPOT-4 (Systeme Probatoire d’Observation de la Terre) at a resolution of 1 km. They found that at the regional scale, the NDVI increased significantly in the growth season, but observed three distinct periods with opposing trends during the entire study period: the NDVI significantly increased by 0.0019 yr⁻¹ in 1982–1994, then decreased by 0.0058 yr⁻¹ in 1994–2000, and again increased by 0.0049 yr⁻¹ in 2000–2011 [15].

Climate drivers and human activities are the main causes of changes in vegetation, and previous studies have examined these factors. Precipitation, temperature, and radiation are the leading climatic drivers that act as the limiting factor for 52%, 31%, and 5% of global vegetated areas, respectively [1]. Climate change has played a major role in the trend of vegetation growth in the Koshi River Basin, and decreases in vegetation were significantly and negatively correlated with a sudden increase in precipitation [15]. The spatial patterns of vegetation growth were varyingly and significantly affected by temperature and precipitation in different periods. Different vegetation types in different regions (such as different climatic regions, different ecological regions, and different soil regions) also have different responses to climate [3]. In Africa and other arid areas, the correlations on NDVI data resulted in the higher correlation between NDVI and precipitation than that with temperature [16]. Forest in Inner Mongolia was positively linked with potential evapotranspiration and temperature, but negatively correlated with precipitation [3]. However, warming-related NDVI at the beginning of the growing season was the main contributor to the NDVI increase during the growing season for forest, meadow steppe, and typical steppe [3]. For semi-arid areas in Iran, yearly precipitation by all vegetation types has close correlations with NDVI fluctuations, rich ranges affected by precipitation, followed by moderate ranges, dry farming and irrigated farming [17]. Vegetation greening in the Himalayas was mainly attributed to human activities such as the use of chemical fertilizers and improved irrigation. In areas where there is no direct human impact, this can be attributed to climate change [14]. The browning in the low-altitude and middle-altitude areas can be attributed to human activities because of the large human disturbance there, whereas the browning in high-altitude areas is most likely caused by climatic factors [14]. The changes in vegetation driven by these two factors also have clear spatio-temporal differences.

Due to the complexity of the relationship between vegetation growth, climate change, and human activities, it is important to understand the status and changes in vegetation as well as reasons for them in order to predict the interactions between terrestrial ecosystems and climate systems at the local scale. In the past, most research on this issue at a regional scale has been based on spatial resolution-related data at the kilometer level—this is not accurate for monitoring the status and change of vegetation at regional and watershed scales scale. In addition, most of previous studies on the influence of climatic factors on vegetation change were based on average temperature, average precipitation, multi-year temperature change, and multi-year temperature precipitation, which does not accurately reflect changes in vegetation in response to climate change compared with the linear time trend method. This study investigates where, when, and why changes in vegetation have occurred in the KRB in the last 18 years using the MODIS (moderate-resolution imaging spectroradiometer) NDVI dataset. This study differs from previous studies in the area due to the use of the high-resolution NDVI dataset (with a resolution of 250 m). We first analyze the status and changes in vegetation in the KRB, and then examine the relationship between the NDVI_GS, and local climate and human activities. Thus, the main influencing factors of vegetation change are obtained. The results can improve our understanding of the underlying drivers of recent changes in vegetation activity in the KRB based on observational analyses, and can help plan for the sustainable development of the local vegetation.

2. Materials and Methods

2.1. Study Area

The KRB is an important transboundary basin across three countries (China, Nepal, and India), covering an area of 87,460.5 km² in the central Himalayas, with geographical coordinates 85°01′ E–88°57′ E and 25°20′ N–29°12′ N (Figure 1). The KRB stretches from the south of the Brahmaputra River in the north to the Nepalese southern frontier in the south, from Kathmandu in the west to the border with Nepal in the east. Its largest drop in elevation, 30 m, is at the junction of the Ganges River with Mount Qomolangma, at an elevation of 8844 m. Records from weather stations in the KRB show that the annual average temperature is 5.1 °C, with an average annual rainfall of 397.6 mm in its Chinese part [18]. The average temperature is 10.5 °C, with an average annual precipitation of 1794.6 mm in the Nepalese part, and the maximum and minimum monthly mean temperatures are 41.8 °C and 15.4 °C, respectively, with an average annual precipitation of 1900 mm in the Indian part. Owing to its high elevation and changes in climate, vegetation across the KRB exhibits clear vertical zonation [15]. The KRB is also one of the world’s most biodiverse areas, rich in flora and fauna, and featuring 6244 species of vascular plants [19].

2.2. Data and Sources

The data used in this study included MODIS NDVI data, meteorological records, and land cover and DEM (digital elevation database) data. The MOD13Q1 NDVI data were downloaded from the National Aeronautics and Space Administration’s website (http://daac.gsfc.nasa.gov/), spanning the period 2000–2018, at a spatial resolution of 250 m and a temporal resolution of 16 days. All images downloaded were mosaicked and re-projected to Albers using the MODIS Reprojection Tools (MRT). To reduce noise induced by poor atmospheric conditions and cloud contamination, we used the Savitzky–Golay filter provided by TIMESAT 3.2 to reconstruct a high-quality NDVI. The meteorological data included data from stations for the period 2000–2018 (five stations in the Chinese KRB), CRU TS 4.03 data (https://www.cru.uea.ac.uk/data), and other interpolated meteorological data, including average temperature, maximum temperature, minimum temperature, and average precipitation. The land cover data for the KRB were obtained from the Land Change and Regional Adaptation Research Group of the Tibetan Plateau, Institute of Geographic Sciences and Natural Resources Research, CAS for 2000–2015 at a spatial resolution of 30 m. The land cover data were reclassified into six classes, i.e., croplands, forests, shrublands, grasslands, sparse vegetation lands, and wetlands. The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global DEM (GDEM) (http://glcfapp.glcf.umd.edu/data/aster/) was used with a spatial resolution of 30 m. To unify the resolution, the gridded meteorological, land cover datasets, and DEM data were re-projected and resampled to a resolution of 250 m.

2.3. Methods

High-quality NDVI data were obtained through preprocessing and Savitzky–Golay filtering [20,21,22]. Based on an analysis of the literature and meteorological data of the KRB, we defined May to October as the growing season each year, and used the R software package to extract the NDVI for this season for each year. The linear time trend method was then used to determine the overall situation of changes in NDVI_GS in the basin. Finally, correlation analysis was used to analyze the trends of climatic factors and the NDVI.

2.3.1. Analysis of Linear Time Trend and Turning Points of MODIS NDVI

A time series analysis of the NDVI was conducted using a linear model (LM) based on the generated NDVI images of the annual growing season from 2000 to 2018 to calculate its linear regression on a spatial grid scale. The trends were quantified by the slope of the regression line. We calculated the significance of each grid. The slope (a) of the regression indicates the mean temporal change, a > 0 represents a positive slope indicating a trend of increase, and a < 0 represents a negative slope, indicating one of decrease [20,21]. The layers of significance (p-values) were categories of the types of change: significant increase (a > 0, p < 0.05), insignificant increase (a > 0, p > 0.05), significant decrease (a < 0, p < 0.05), and insignificant decrease (a < 0, p > 0.05) [20,21].

The BFAST (Breaks for Additive Seasonal and Trend) model was used to analyze the turning points, which was proposed by Verbesselt in 2010 [23]. It was initially used to identify vegetation disturbance using remote sensing data. The algorithm is widely applicable and can be used in meteorological, hydrological, economic, and other fields [24,25]. The long-term trend T_i was piecewise fitted using a linear model, for each segment, t_i < t ≤ t_i+1, t₀ = 0 is defined. The fitted linear model is as follows:

Tt = αi + βit,

(1)

where i is the location of the mutation point, with the value of 1 to p, there are a total of p mutation points. The intercept α_i and slope β_i of the mutation point and the linear mode on both sides can analyze the degree and direction of the mutation and the gradient slope [25].

The algorithm of BFAST is implemented in the BFAST software package of R language.

2.3.2. Correlation between NDVI and Climate Variables

Because the linear model can show only independent contributions and superimposed relationships between the variables, to quantify the influence of individual climate drivers (precipitation and temperature), we examined the correlations between NDVI_GS and them. The Pearson correlation coefficient was used in our model to reflect the effect of interactions between the total annual precipitation trend (PRE), annual average temperature trend (TMP), annual maximum temperature trend (TMX), annual minimum temperature trend (TMN), and NDVI_GS trend in 2000–2018 at the pixel scale. The value of Pearson correlation was in the range (−1, 1). The larger the value was, the stronger the relationship was between vegetation and the driving factors [22]. The formula for this is:

$${\mathrm{r}}_{\mathrm{xy}}=\frac{{{\displaystyle \sum }}_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2\times {{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}},$$

(2)

where r_xy is the correlation coefficient between NDVI_GS and climatic variables, x_i represents the NDVI_GS value in the ith year, y_i represents the climatic variables in the ith year, and $\overline{x}$ and $\overline{y}$ represent the averages of NDVI_GS and the climatic variables, respectively [22].

3. Results

3.1. Spatio-Temporal Change in NDVI_GS from 2000–2018 and the Turning Points

The average NDVI_GS for 2000–2018 is shown in Figure 2a. Areas featuring glaciers and bodies of water have been left blank. The range of average NDVI_GS in 2000–2018 was from 0.1 to 0.9. NDVI_GS high value (0.5–0.8) areas with dense vegetation cover were found over the Nepalese and Indian parts of the KRB, accounting for 58% of all grid cells. These areas mainly featured forests and cropland. NDVI_GS low value (0.1–0.2) areas reflecting poor vegetation cover were observed over the Chinese part of the KRB, accounting for 18% of all grid cells. These areas were mainly composed of grassland cover. The area with NDVI_GS value from 0.9–1 was the smallest. The second smallest area was NDVI_GS value from 0.3–0.4.

The inter-annual variation (average over all grids cells using the mean NDVI_GS) of NDVI_GS from 2000–2018 is shown in Figure 2b. The average NDVI_GS showed a significant (p < 0.001) upward trend with a greening rate of 0.0023 yr⁻¹ in 2000–2018. The annual NDVI_GS value increased by 9.8% over the 18-year period. The highest annual NDVI_GS was 0.52 in 2017, 0.03 higher than the multi-year average NDVI_GS (0.49), followed by that in 2018 (0.51). The lowest annual value, 0.46, occurred in 2000 and 2004.

The NDVI_GS increased across 78% of the KRB in 2000–2018 (83% of the grid cells) (Figure 2c). Overall, 27% of the grid cells showed trends of significant (significance, p < 0.05) increase, particularly in the Indian KRB, which mainly contained cropland. It also increased in the valleys in Nepal, covered by low-altitude cropland and forests. Only 1.1% of the grid cells showed trends of significant (p < 0.05) decrease, especially around Kathmandu, the capital of Nepal, where cropland had been converted into construction land. Significant decreases were also observed in the Chinese KRB, at the juncture of Dinggye County and Gamba County. Insignificant increases accounted for 56%, and were mainly distributed in the Chinese and Nepalese KRB. Insignificant decreases (accounting for 16%) mainly occurred in northern part of the Nepal KRB and the northeastern part of the Chinese KRB.

At the regional scale, the TPs (Turning points) of NDVI_GS in 2000–2018 in the MODIS NDVI dataset occurred in 2004 and 2007 (Figure 2d). Although a statistically significant trend was observed for NDVI_GS over the entire study period, there were three distinct periods with different trends of variation: 2000–2004, 2004–2007, and 2007–2018 (Figure 2d). NDVI_GS showed an increasing trend in the three periods with different annual growth rates. In 2000–2004, NDVI_GS spatially averaged over the KRB marginally increased by 0.0025 yr⁻¹, which was close to the overall average annual growth rate. In 2004, NDVI_GS suddenly decreased by 0.03, and showed a rapid growth at an annual rate of 0.004 yr⁻¹ in 2004–2007, which was the fastest growing period of NDVI_GS in KRB. In the third period, NDVI_GS again suddenly decreased by 0.025 in 2007, and then increased by 0.0011 yr⁻¹ in 2007–2018, which was the slowest growing period of NDVI_GS. Therefore, although the NDVI_GS in KRB generally showed an increasing trend, its greening rate was slowing down, especially in the recent ten years.

3.2. Variations in NDVI_GS along Elevation Gradient

The average NDVI_GS value generally decreased with elevation (Figure 3a). In more detail, it can be divided into three distinct stages. The first stage is the high-value sector below 3600 m, where the average NDVI_GS was greater than 0.7. Because of the particularity of the huge altitude difference in the study area, a large area of croplands is distributed in the plain area less than 100 m, and the high NDVI_GS value of croplands leads to the average value NDVI_GS of low altitude areas being higher. The second stage is the mid-value sector from 3600 m to 4100 m, where the average NDVI_GS was 0.48. The third stage was the low value sector above 4100 m, with an average NDVI_GS of 0.12.

The annual average change rate (ACR) of NDVI_GS was positive overall, which showed that each elevation gradient presented an upward trend (Figure 3a). However, the ACR value of NDVI_GS generally decreased with elevation: it decreased from 0.0064 yr⁻¹ at 30 m to 0.0019 yr⁻¹ at 1400 m, gradually increased to 0.0055 yr⁻¹ at 2300 m, sharply decreased to 0.0011 yr⁻¹ at 3400 m, sharply increased to 0.0036 yr⁻¹ at 3900 m, and then decreased rapidly to 0.0013 yr⁻¹ at 4500 m. The curve of the ACR showed two prominent peaks and two troughs. One peak of 0.64% occurred at 0–100 m and the other of 0.0055 yr⁻¹ at 2200–2300 m. The first trough was 0.0019 yr⁻¹ at 1300–1400 m and the second was 0.0011 yr⁻¹ at 3300–3400 m.

The ACR value was divided into greening and browning rates for further analysis (Figure 3b). Large greening areas distributed in the ranges of elevation of 30 m to 1800 m and 4300 m to 5500 m, particularly at elevation lower than 100 m covered with croplands, where the greening area accounted for 46% of the total. The greening rate decreased with elevation below 900 m (from 0.0068 yr⁻¹ to 0.0039 yr⁻¹), gradually increased with the elevation from 900 m to 3400 m (0.0078% yr⁻¹), and then rapidly decreased with elevation from 3400 m to 4300 m (0.0024 yr⁻¹); above 4300 m, the greening rate was relatively stable. Two peaks were observed in greening rate, one was 0.0068 yr⁻¹ at 0–100 m and another was 0.0078 yr⁻¹ at 3300–3400 m. The browning rate decreased with elevation below 1000 m (from −0.014 yr⁻¹ to −0.0038 yr⁻¹), gradually increased with elevation from 1000 m to 3100 m (−0.012 yr⁻¹), and then decreased from 3100 m to 4500 m (−0.0017 yr⁻¹); above 4500 m, the browning rate changed smoothly. Two troughs were observed, one was −0.014 yr⁻¹ at 0–100 m and another one was −0.011 yr⁻¹ at 2900–3000 m.

3.3. Variation in NDVIGS with Land Cover Type

Through the integration and reclassification of land cover data in KRB, vegetation in it was divided into six types: cropland, forest, shrubland, grassland, sparse vegetation, and wetland. The highest average NDVI_GS was obtained for forests (0.76), followed by shrublands (0.64) and croplands (0.61) (Figure 4a), but in terms of vegetation area, the area of croplands was the largest, followed by forests and grasslands. From the perspective of the ACR of NDVI_GS (Figure 4a), the annual NDVI_GS also exhibited trends of increasing from 2000 to 2018 in different vegetation types, with greening rates ranging from 0.0014 yr⁻¹ to 0.0059 yr⁻¹. The increase rate of croplands was the fastest, followed by forests, shrublands, and wetlands.

Croplands covered the largest greening area (approximately 20,423.33 km²), followed by forests (5525.27 km²) and grasslands (3667.16 km²) (Figure 4b)—those constituted the main force of NDVI greening in KRB. In terms of the change rate, the fastest greening rate was also observed for croplands (0.0062 yr⁻¹), followed by shrubland (0.0048 yr⁻¹) and forests (0.0043 yr⁻¹) (Figure 4b). The NDVI_GS increased weakly in the natural ecosystems, but strongly in cropland [3], which is the result of human land use activities. The vegetation type with largest browning area was grasslands (approximately 578.65 km²), which had the lowest browning rate (−0.0022 yr⁻¹) (Figure 4b). The second largest browning area was in the croplands (465.23 km²), which had the highest browning rate (−0.0076 yr⁻¹), identical to that of shrublands (Figure 4b). The area of KRB greening was far larger than that of browning, leading to an overall trend of greening.

The variations in NDVI_GS at different elevation zones were highly correlated with vegetation type (Figure 4c,d). The annual NDVI_GS values for all vegetation types showed relatively stable, and fluctuation and a slight increase occurred when elevation was below 3000 m, with grasslands recording the highest fluctuation. When elevation was above 3000 m, the annual NDVI_GS of all vegetation types decreased due to the rapid decrease of vegetation area with high NDVI_GS value. Wetlands, grasslands, and forests in particular decreased sharply with elevation above 3000 m, whereas the NDVI_GS of croplands showed a slight rise at about 4500 m, that is, the distribution of cropland like highland barley and oat grass in high elevation areas of Tibet. Sparse vegetation distributed in the high elevation area above 4000 m, and the NDVI_GS of sparse vegetation began to decrease sharply at about 4400 m. In terms of change rates, croplands, forests, and shrublands were relatively stable, showing an increase in greening rate with elevation. However, when above 2500 m, greening rates of them began to decrease. For instance, the change rate of grasslands changed considerably with elevation—less than 1000 m showed a greening trend, between 1100 m to 1500 m showed a browning trend, and between 1500 m to 3100 m showed a greening trend, then between 3100 m to 3400 m showed a browning trend, and above 3400 m showed a greening trend.

3.4. Correlation between Trends of NDVI_GS and Climate Factors

Correlation between Trend of NDVI_GS and Climate Change

There were clear regional differences in terms of correlation between changes in NDVI_GS and climatic factors (Figure 5). Regions with significant correlations between changes in NDVI_GS and climate change accounted for 38.7% of the total area of the KRB. However, there was little overlap between the regions of significant correlation and the regions of significant increase and decrease in NDVI_GS. Therefore, climate change may not be the dominant factor in vegetation change in the region. There was an insignificant positive correlation between the PRE and NDVI_GS in the Chinese part of KRB that occupies a large area where NDVI_GS showed an insignificant increase trend. This indicates that there was an insignificant positive correlation between vegetation increase and PRE in Chinese part of KRB. Areas with a significant positive correlation between changes in PRE and NDVI_GS occupied a small area, and were mainly distributed in grassland vegetation around rivers, lakes, and wetlands in Nyalam county, Dingri county, Sakya County, and Dinggye county of China, where NDVI_GS also showed insignificant increase trend. An insignificant positive correlation was observed between the TMP and grasslands in the west of Chinese part of KRB, and an insignificant negative correlation between the TMP and grassland in the middle-east of Chinese part of KRB.

An insignificant negative correlation was observed between vegetation change and the PRE in Nepalese part of KRB, where NDVI_GS showed an insignificant decrease and mainly covered high altitude croplands and forests. The significant increase of NDVI_GS had an insignificant positive correlation with PRE in low- and medium-altitude areas covered with forests and croplands. In addition, at low altitudes, forest can be more sensitive to variations in precipitation because of the larger growing season. Changes in NDVI_GS and the TMP exhibited an insignificant negative correlation, and there was an insignificant positive correlation between the TMP and forests at high elevation. In addition, at high altitudes, the reduction in evapotranspiration makes vegetation more sensitive to temperature and the length of the growing season. An insignificant negative correlation was observed between the TMX and forests and croplands in the Nepalese part, and an insignificant negative correlation between changes in NDVI_GS and the TMN. The smaller the TMN was, the larger the changes in NDVI_GS were. Most of the croplands and forests were negatively correlated with the TMN.

No significant positive correlation was observed between cropland and the PRE in the Indian part of the KRB. Most of the cropland there had a negative correlation with the TMP that was not significant. NDVI_GS and TMP were positively correlated in the flood area on the right side of the river, but the correlation was not significant. There was no significant negative correlation between croplands and the TMX in India, and a significant negative correlation between local areas (on the west side of the border between Nepal and India). No significant positive correlation was observed between areas on either sides of the river and the TMX, and a significant negative correlation was obtained in the Indian KRB between the TMN and changes in the NDVI. Most of the regions were negatively correlated.

A significant positive correlation was obtained between the NDVI_GS and the TMX around the glaciers along the border between Nepal and China. Most of these areas are glacial retreat areas, which are positively correlated with the maximum temperature, and with the increase of the maximum temperature, the glacier retreat rate is faster.

4. Discussion

4.1. Trends of Temporal and Spatial Changes in NDVI_GS

In this study, we investigated the spatio-temporal variations in NDVI_GS over the KRB in 2000–2018 through trend analysis using MODIS NDVI, and examined their correlations with climatic factors. The results show a greening trend (greening rate was 0.0023 yr⁻¹) from 2000–2018 over the KRB, which is consistent with the previous study by Zhao et al. [1], where the global mean NDVI_GS showed an increasing trend from 2000–2013, with significant greening over Eurasia. Zhang et al. [15] also found a significant increase (0.0034 yr⁻¹) in the average NDVI_GS in the KRB in China and Nepal in 2000–2011, but the greening rate was higher than this study. The NDVI data resolution may be one reason—MODIS NDVI data with a resolution of 250 m was used in this study, while Zhang et al. used SPOT-4 Vegetation (VGT) data with a resolution of 1 km. The second reason can be explained by the TPs analysis of this study. The TPs of NDVI_GS in 2004 and 2007 in the KRB were consistent with the result obtained by Ding et al. [26], which showed that a transition in the phenological index occurred at the end of the growing season in 2004–2007 [26]. Although the overall NDVI_GS of the KRB showed an increasing trend, the greening rate of NDVI_GS slowed down after the second TPs with a greening rate of 0.0011 yr⁻¹ from 2007–2018 lower than the overall increasing rate.

In the context of trends of spatial changes in NDVI_GS, our results show a greening trend from 2000 to 2018 over the Indian KRB, mainly covered by cropland, and at low elevation in Nepal, featuring croplands and forests. This was consistent with the results of many previous studies [12,14,27,28,29]. For instance, de Jong et al. [12] found that parts of the Indian KRB exhibited a greening trend for 20 years or longer (1982–2008). Wang et al. [27] also found widespread vegetation-related greening in India. Chen et al. [30] found that the KRB in China and India led in terms of global greening, with the greening in China mainly due to forests (42%) and croplands (32%), and in India mostly from croplands (82%) with a minor contribution by forests (4.4%). This is consistent with the results of this study. Mishra and Mainali [14] showed that a majority of Himalayan vegetation (include parts of Nepal) has been undergoing greening since 2000. The area of NDVI_GS trend toward browning in KRB occupied a small area, mainly occurring around Kathmandu, where croplands had been rapidly transformed into construction land from 2000 onward. Significant areas of decrease were also found in the Chinese part of the KRB at the juncture of Dinggye county, which was consistent with results by Han et al. [28]. Based on NDVI data, Liu et al. [29] studied the grassland changes in the Tibetan Plateau, and also found that the grassland with a downward trend was mainly distributed in Tibet region, while the grasslands that experienced upward trends were mainly in the Qinghai Province.

4.2. Response in Terms of NDVI_GS to Climate Change and Human Activities

Changes in vegetation are mainly determined by climate change [31,32,33] and human activities [31,34,35,36,37]. Many studies have shown that precipitation is the main factor influencing an increase in NDVI_GS [38,39,40,41]. This is not consistent with the results of this study. An insignificant positive correlation was noted between the PRE and NDVI_GS over a large area in the Chinese KRB, whereas the changes in NDVI_GS there showed insignificant increases. A large part of the area was covered with grasslands owing to local precipitation [15,42,43] and protective governmental policies, such as ecological protection projects [44,45] and nature reserves [46,47,48], and this led to an increase in NDVI_GS. The research of Zhang et al. [15] in China and India part of KRB pointed out that climate change may be the main cause of vegetation growth in this area, the decrease of NDVI in the periods of 1982–1994 and 1994–2000 was significantly and negatively correlated with the sudden increase in precipitation. However, the study period for vegetation change being significantly correlated with precipitation was different from the research period of this study. Thus, the results are not comparable due to the spatial and temporal differences in the impacts of climate change and human activities on vegetation. Second, the above literature focuses on the impact of climate change on vegetation, without considering the impact of human activities, especially in the context of frequent human activities in recent years.

The NDVI_GS of croplands and forests in the Indian and Nepalese KRB increased significantly, those also were not significantly related to climatic factors. Although the increases in carbon dioxide concentration and nitrogen deposition in the atmosphere were the most likely climatic causes of greening [14], the success of sustainable forestry practices (community forestry in Nepal) [49,50], and the increase in agricultural fertilization [51,52] and irrigation facilities [53,54,55,56] in Nepal and India might have been the driving force for it [14]. A rare instance of the success of community forestry was noted in Nepal, especially at low elevation (200–2000 m) [57]. Wang et al. [27] found that enhanced moisture conditions over dry regions, coinciding with a decline in monsoon, were mainly responsible for the widespread trend of greening. Chen et al. [30] found that food production in India has increased since 2000 mostly due to an increase in the harvested area through multiple crops facilitated by fertilizer use and surface and groundwater irrigation. Water supply was the main factor restricting vegetation growth in dry seasons and areas. However, the spatial and temporal distributions of precipitation in India were uneven. High-intensity precipitation was confined in time to the rainy season (July to September) when the monsoon comes, which most likely led to floods that directly destroyed vegetation and crops. As a result, the vegetation coverage in the flood area is reduced instantaneously [58], while after the withdrawal of the southwest monsoon, there is almost no precipitation in the dry season and the growth of vegetation is restricted. The uneven spatial and temporal distribution of precipitation is also an important factor limiting vegetation growth.

The significant decrease in NDVI_GS in the KRB was mainly distributed around Kathmandu in Nepal and Gangba county in China. Owing to the rapid urbanization in Kathmandu from 2000 to 2018 [59,60], most cropland around the city has been converted into construction land [60], which led to a decline in the NDVI around it [61,62]. There was no significant relationship between NDVI_GS and climatic factors in Gangba county. As a typical animal product of Tibet, Gamba sheep obtained the registration certificate of geographical indication of agricultural products of the People’s Republic of China in December 2014, and registered the trademark of Gamba sheep in 2015, which further improved the brand awareness and market price of Gamba sheep [63,64]. In addition, Gangba county has an extremely fragile ecological environment, and grazing pressure has led to a degradation of the grassland [65,66]. Liu et al. [29] also found the decreasing trend of grassland in Tibet region mainly related to the increasing number of Tibetan people and livestock. Therefore, the government should attend to measures for its ecological protection.

Changes in NDVI_GS around glaciers at the junction of China and Nepal had a significant positive correlation with the TMX. Perhaps the supply of glacial meltwater promoted the growth of vegetation, but NDVI_GS decreased insignificantly farther from glaciers. The reverse showed the contribution made by glacial meltwater supplements. Anderson et al. [67] studied the expansion in vegetation in the Himalayas, and concluded that the strongest and most significant trends were at elevations of 5000–5500 m across the Hindu Kush and the Himalayas, covering the junction of China and Nepal. At higher elevations, flatter areas exhibited stronger trends [67].

5. Conclusions

This study showed that vegetation has significantly increased over the KRB from 2000 to 2018, which is consistent with previous results noted across Asia, especially in Nepal and India [15,62]. The NDVI_GS of cropland in the Indian part of the KRB, and of cropland and forest in its Nepalese part showed a trend of significant increase. Climate change did not play a major role in an increase in NDVI_GS in the KRB. Changes in NDVI_GS were mainly affected by human activities and climate change, and the influence of the former in low-altitude areas was more significant. Increases in NDVI_GS in croplands in India and Nepal occurred mainly owing to human interference, such as improvements in agricultural fertilization and irrigation. The success of community forestry in Nepal had a positive effect on the increase in the NDVI_GS of local forests. Climate also had a certain impact on the increase in NDVI_GS, such as a significant negative correlation between it and the TMN in the Indian KRB, and insignificant positive correlations with the PRE. The decrease in NDVI_GS around Kathmandu was mainly due to urbanization. Increases in NDVI_GS in the KRB were thus mainly affected by human activities, and by climate change to some extent.