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

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 NDVIGS (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 NDVIGS of the local forest. Climate also had a certain impact on the increase in NDVIGS. A significant negative correlation was observed between NDVIGS 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). NDVIGS significantly decreased over a small area, mainly around Kathmandu, due to urbanization. Increases in NDVIGS in the KRB have thus been mainly affected by human activities, and climate change has helped increase it to a certain extent.


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

Study Area
The KRB is an important transboundary basin across three countries (China, Nepal, and India), covering an area of 87,460.5 km 2 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].
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 15 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 NDVIGS, 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.

Study Area
The KRB is an important transboundary basin across three countries (China, Nepal, and India), covering an area of 87,460.5 km 2 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].

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.

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.

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 = 0 is defined. The fitted linear model is as follows: 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.

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: 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 x and y represent the averages of NDVI GS and the climatic variables, respectively [22].

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 −1 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)  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.
Sustainability 2020, 12, x FOR PEER REVIEW 6 of 15 which was the slowest growing period of NDVIGS. Therefore, although the NDVIGS in KRB generally showed an increasing trend, its greening rate was slowing down, especially in the recent ten years.  Figure 2b shows coefficients of the slope, and (p) is the significance test.

Variations in NDVIGS along Elevation Gradient
The average NDVIGS 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 NDVIGS 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 NDVIGS value of croplands leads to the average value NDVIGS of low altitude areas being higher. The second stage is the mid-value sector from 3600 m to 4100 m, where the average NDVIGS was 0.48. The third stage was the low value sector above 4100 m, with an average NDVIGS of 0.12.
The annual average change rate (ACR) of NDVIGS was positive overall, which showed that each elevation gradient presented an upward trend (Figure 3a). However, the ACR value of NDVIGS generally decreased with elevation: it decreased from 0.0064 yr −1 at 30 m to 0.0019 yr −1 at 1400 m, gradually increased to 0.0055 yr −1 at 2300 m, sharply decreased to 0.0011 yr −1 at 3400 m, sharply increased to 0.0036 yr −1 at 3900 m, and then decreased rapidly to 0.0013 yr −1 at 4500 m. The curve of  Figure 2b shows coefficients of the slope, and (p) is the significance test.

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. rate was relatively stable. Two peaks were observed in greening rate, one was 0.0068 yr −1 at 0-100 m and another was 0.0078 yr −1 at 3300-3400 m. The browning rate decreased with elevation below 1000 m (from −0.014 yr −1 to −0.0038 yr −1 ), gradually increased with elevation from 1000 m to 3100 m (−0.012 yr −1 ), and then decreased from 3100 m to 4500 m (−0.0017 yr −1 ); above 4500 m, the browning rate changed smoothly. Two troughs were observed, one was −0.014 yr −1 at 0-100 m and another one was −0.011 yr −1 at 2900-3000 m.

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 NDVIGS 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 NDVIGS (Figure 4a), the annual NDVIGS also exhibited trends of increasing from 2000 to 2018 in different vegetation types, with greening rates 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 −1 at 30 m to 0.0019 yr −1 at 1400 m, gradually increased to 0.0055 yr −1 at 2300 m, sharply decreased to 0.0011 yr −1 at 3400 m, sharply increased to 0.0036 yr −1 at 3900 m, and then decreased rapidly to 0.0013 yr −1 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 −1 at 2200-2300 m. The first trough was 0.0019 yr −1 at 1300-1400 m and the second was 0.0011 yr −1 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 −1 to 0.0039 yr −1 ), gradually increased with the elevation from 900 m to 3400 m (0.0078% yr −1 ), and then rapidly decreased with elevation from 3400 m to 4300 m (0.0024 yr −1 ); above 4300 m, the greening rate was relatively stable. Two peaks were observed in greening rate, one was 0.0068 yr −1 at 0-100 m and another was 0.0078 yr −1 at 3300-3400 m. The browning rate decreased with elevation below 1000 m (from −0.014 yr −1 to −0.0038 yr −1 ), gradually increased with elevation from 1000 m to 3100 m (−0.012 yr −1 ), and then decreased from 3100 m to 4500 m (−0.0017 yr −1 ); above 4500 m, the browning rate changed smoothly. Two troughs were observed, one was −0.014 yr −1 at 0-100 m and another one was −0.011 yr −1 at 2900-3000 m.

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 −1 to 0.0059 yr −1 . The increase rate of croplands was the fastest, followed by forests, shrublands, and wetlands.
browning, leading to an overall trend of greening.
The variations in NDVIGS at different elevation zones were highly correlated with vegetation type (Figure 4c,d). The annual NDVIGS 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 NDVIGS of all vegetation types decreased due to the rapid decrease of vegetation area with high NDVIGS value. Wetlands, grasslands, and forests in particular decreased sharply with elevation above 3000 m, whereas the NDVIGS 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 NDVIGS 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.  Croplands covered the largest greening area (approximately 20,423.33 km 2 ), followed by forests (5525.27 km 2 ) and grasslands (3667.16 km 2 ) (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 −1 ), followed by shrubland (0.0048 yr −1 ) and forests (0.0043 yr −1 ) (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 2 ), which had the lowest browning rate (−0.0022 yr −1 ) (Figure 4b). The second largest browning area was in the croplands (465.23 km 2 ), which had the highest browning rate (−0.0076 yr −1 ), 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.

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.

Trends of Temporal and Spatial Changes in NDVIGS
In this study, we investigated the spatio-temporal variations in NDVIGS over the KRB in 2000-2018 through trend analysis using MODIS NDVI, and examined their correlations with climatic 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.

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

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

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.
Author Contributions: Y.Z. and X.W. had the original idea and designed the research. X.W. performed the study, processed and analyzed the data and wrote the manuscript. X.S., Z.W. and B.P. had insights on the revision of the manuscript and suggestions for improvement. Y.Z. came up with the theme and assisted with the suggestion of the study area, provided project and financial support, and had insights on the revision of the manuscript and suggestions for improvement. Q.L. and B.Z. provided help and guidance for data processing and model adjustment. F.X. interpreted the land cover data of KRB in 2015. All authors participated in writing the final version of the manuscript, and reviewed and approved it. All authors have read and agreed to the published version of the manuscript.