Assessing Climate Change Trends and Their Relationships with Alpine Vegetation and Surface Water Dynamics in the Everest Region, Nepal

The Himalayas, especially the Everest region, are highly sensitive to climate change. Although there are research works on this region related to cryospheric work, the ecological understandings of the alpine zone and climate impacts are limited. This study aimed to assess the changes in surface water including glacier lake and streamflow and the spatial and temporal changes in alpine vegetation and examine their relationships with climatic factors (temperature and precipitation) during 1995–2019 in the Everest region and the Dudh Koshi river basin. In this study, Landsat time-series data, European Commission’s Joint Research Center (JRC) surface water data, ECMWF Reanalysis 5th Generation (ERA5) reanalysis temperature data, and meteorological station data were used. It was found that the glacial lake area and volume are expanding at the rates of 0.0676 and 0.0198 km3/year, respectively; the average annual streamflow is decreasing at the rate of 2.73 m3/s/year. Similarly, the alpine vegetation greening as indicated by normalized difference vegetation index (NDVI) is increasing at the rate of 0.00352 units/year. On the other hand, the annual mean temperature shows an increasing trend of 0.0329 °C/year, and the annual precipitation also shows a significant negative monotonic trend. It was also found that annual NDVI is significantly correlated with annual temperature. Likewise, the glacial lake area expansion is strongly correlated with annual minimum temperature and annual precipitation. Overall, we found a significant alteration in the alpine ecosystem of the Everest region that could impact on the water–energy–food nexus of the Dudh Koshi river basin.


Introduction
Climate is a major driving factor for the distribution of the global forest ecosystem [1]. Climate, especially temperature and precipitation, varies along latitude and continental position [2]. Based on temperature and precipitation, the Köppen-Trewartha classification has divided global climate into six groups: tropical, dry climate, sub-tropical, temperate, boreal, and polar climate [3]. Besides this, altitude also creates a mountain system which is defined by a minimum elevation of 300 m and minimum slope of 2 • over 25 km [4]. The mountain system contains distinct thermal belts, which are distinguished into a mountain (>15 • C), montane (6.4-10 • C), alpine (3.5-6.4 • C), and Nival belts (<3.5 • C) [5]. In the broader context, the mountain system can be grouped into three belts: montane belt, alpine belt, and Nival belt. Likewise, the mountain system of Nepal is also divided into six life Besides vegetation, climate warming has also been found to affect cryospheric processes (such as changes in snow cover and glacier permafrost) and related hydrological processes (such as water cycle and balance and river discharge) [7]. For example, globally, it was reported that during 1990-2018 the number and area of glacial lakes have rapidly increased from 9414 (5930 km 2 ) to 14,394 (8950 km 2 ), and the glacial lake volume has increased from 105.7 to 156.5 km 3 [24]. Similarly, the increasing number and area of glacier lakes were also revealed in the central Himalayas region during 1990-2010 [25]. At the scale of Nepal Himalaya, the number and area of glacier lakes have increased from 606 (55.53 km 2 ) to 1541 (80.95 km 2 ) from 1977 to 2017 [26]. It was revealed that the glacier has retreated at the rate of 10-59 m/year in the Everest region [27]. The changes in the glacial lakes will favor vegetation growth by increasing the soil moisture and temperature. However, the glacial lakes at a higher altitude create a threat of glacial lake outburst flood (GLOF) to lowland people. Hence, the trend of glacial lake expansion by area and volume should be monitored on a long-term basis.
Although there are some studies on the vegetation dynamics Koshi river basin, and over Nepal, the study of high altitude forest ecosystem in response to climate change is still unexplored. This region has less human disturbance for land-use change because of extreme climate. Because of this, the alpine ecosystem is considered a quick indicator of the impacts of climate change. The long-term analysis of alpine vegetation of the Everest region will provide a better understanding of climate change interaction with high altitude. In addition, the trend analysis of the number and area of glacial lakes over a long-term basis is also new and useful for the prediction of water stress, adaptability, vulnerability, and overall water resource management in Nepal Himalaya. Specifically, to understand the relationships between climate variables and their potential impacts on the Everest region of Nepal, this study aimed to answer three research questions: (i) What are the trends of temperature and precipitation in the Everest region? (ii) What are the trends of glacial lake and streamflow changes in the Everest region? (iii) What are the NDVI trends and areal changes of alpine vegetation in the Everest region?

Site Description
This study focused on the high mountain system of Nepal in the Mount Everest region. This region is also designated as Sagarmatha National Park covering an area of 1413 km 2 . According to Garrard et al. (2016), regarding the land use/land cover of Sagarmatha National Park in 1992, bare rock/soil plus permanent ice covers 985.26 km 2 (69.73% of total land area); rangeland, 301.49 km 2 (21.34%); forest, 110.27 km 2 (7.8%); agriculture and settlement, 10.21 km 2 (0.72%); lake, 5.77 km 2 (0.41%) [28]. According to another data source from ICIMOD (The International Center for Integrated Mountain Development, an intergovernmental institution for the study of Hindu Kush Himalaya), the land cover in 1992 of the Sagarmatha National Park comprised snow (590.38 km 2 ), rangeland (301.84 km 2 ), forest (120.29 km 2 ), water (5.70 km 2 ), and other lands (394.79 km 2 , Figure 1) (https://rds.icimod.org (accessed on 30 July 2021)).  Figure 2 shows the elevation zone of the study area covering the Dudh Koshi watershed area. This map was prepared based on the life zone of Nepal using QGIS software. Based on elevation, the study area can be divided into temperate (2845-3000 m), alpine (3000-5000 m), and nival zone (>5000 m). The alpine vegetation generally extends from 3000 to 5000 m in elevation. The streamflow is based on the accumulated flow from whole watershed areas of the Dudh Koshi river basin including the study area ( Figure 2).

Climate Data
The basic technique to understand climate change is to analyze long-term records of climate variables from atmospheric observations. To check the seasonal variability, the annual climate is differentiated into four seasons: winter (December, January, and February; DJF), spring (March, April, and May; MAM), summer (June, July, and August; JJA), and autumn (September, October, and November; SON). However, the period of the season varies with continental position; thus, Nepal has adjusted the seasons into winter (DJF), spring (MAM), summer (JJAS), and autumn (ON), especially for precipitation and hydrological analysis [29].
In this study area, there is one climate station for observation of precipitation in Chaurikharka and one hydrological station on the downstream point of the Dudh Koshi watershed area at Rabuwa Bazar, operated by the Department of Hydrology and Meteorology (DHM), Nepal, as shown in Figure 2 [30]. However, the DHM does not have a temperature station in the vicinity of the Everest region. Alternatively, the temperature was analyzed by using ERA5 reanalysis data. The details of meteorological stations are listed in Table 1. The handling of missing data is also important in the quality control of climate data. The threshold of omission of observation data for analysis seems different among the researchers. Zhai et al. (1999) omitted stations with more than 20 precipitation data missing in a year and 3 consecutive missing years or 5 randomly missing years [31]. Zubieta et al. selected a series with less than 5% missing data in a month [32]. Chapman and Walsh [33], Klein Tank et al. [34], and Haylock et al. [35] set a threshold of no more than 20% of missing daily, monthly, and annual data series, respectively. In this study, the proportions of monthly missing data in precipitation and streamflow were 1.35% and 5.56%, respectively.

Temperature
The ERA5 reanalysis temperature data (1995-2019) were extracted from Climate Data Store, Copernicus (https://cds.climate.copernicus.eu) (accessed date 30 July 2021), which have been produced by ECMWF/Copernicus Climate Change Service using R software. These reanalysis data are available on a regular latitude-longitude grid of 0.25 • × 0.25 • resolution, equivalent to approximately 31 km × 31 km, from 1950 to the present [36]. In this study, the ERA5 temperature data at Khumjung points were interpolated from gridded ERA5 data. ERA5 data are widely used to monitor climate change even inaccessible places and have been found effective and applicable [37,38]. Since reanalysis data are complete and consistent in time, missing data and outliers are not found in time-series data. After acquiring, the homogeneity of data was examined using the RHTests software package, which is based on penalized maximal F test [39]. The monthly data were converted into four seasonal series (winter (DJF), spring (MAM), summer (JJA), and autumn (SON)) and one annual series by the averaging method. The trend analysis of annual and seasonal time-series data was performed by using the modified Mann-Kendall test by the Hamed and Rao approach [40]. This non-parametric Mann-Kendall test has been used to evaluate the monotonic upward or downward trend of variables over time [41]. The Hamed and Rao approach has been used to address the autocorrelation in observed climate data [40]. Additionally, to determine the magnitude of change and forecasting, time-series data were also analyzed by linear regression test.
Reanalysis data are the best possible estimate and evenly spaced data grid of the current state of the Earth system using all available observations with a weather forecasting model. However, it is essential to validate with observation data to ensure their reliability. In this study, the ERA5 reanalysis data were compared with monthly data of 1995 (12 samples, 4% of population) of the nearest observation station, the Pyramid Observatory Laboratory (5050 m), which is located at 27.9592 • N, 86.8131 • E, Khumbu Glacier, in the Sagarmatha National Park, Nepal [42]. The independent t-test was used to compare the significance of the mean difference within two independent data.

Precipitation
Chaurikharka weather is the nearest permanent weather station, managed by DHM, Nepal. DHM is using a traditional US standard (8" diameter) manual rain gauge to measure the rainfall form of precipitation [43]. For a better understanding of climate variation, climate observation stations should be placed at adequate spatial resolution. According to World Meteorological Organization (WMO), it is recommended that the spacing of climatological stations should not be more than 100 km [44]. The study area extends around 20 km east-west and 50 km north-south. Within this criterion, the Chaurikharka station is very likely to represent the study area.
The daily and monthly precipitation data (1995-2019) were obtained from DHM, Nepal. Then, collected data were scanned to detect outlier and missing data. The outlier data were sliced by the interquartile range (IQR) method. An outlier was sliced out when X i − q 50 > f × IQR; where X i is ith observation, q 50 is the 50th percentile, f is the multiplication factor (default f = 3), and IQR = q 75 − q 25 [45]. The proportion of missing data in the monthly series was found to be 1.35%. These missing monthly data were filled by the all-time (1995-2019) mean value of the corresponding month's precipitation data. For example, a missing January datum was filled by the mean value of all January data from 1995 to 2019. After this, the homogeneity of data was examined using the RHTests software package. Thereafter, the monthly data were converted into four seasonal series (winter (DJF), spring (MAM), summer (JJAS), and autumn (ON)) by the summation method. The trend of seasonal and annual time-series data was analyzed by the modified Mann-Kendall test. Additionally, the linear regression test was used to evaluate the trendline of time-series data.

Glacial Lakes
The glacier lake area data (1995-2019) were extracted from JRC Global Surface Water Dataset, developed by The European Commission's Joint Research Centre [46]. This dataset has already been validated using ground control points. The extent and change of water surface are provided at 30 m spatial resolution. The surface area of glacier lakes by month was computed and extracted using Google Earth Engine [47]. The data were further converted into seasonal and annual series by the averaging method. The trend analysis of surface area changes was performed by the modified Mann-Kendall test and linear regression test.
The volume of the glacier lake was determined using the area-volume empirical relationship developed by Huggel et al. [48] with the equation: V = 0.104 × A 1.42 , where V is the volume of glacier lake (m 3 ) and A is the surface area of glacier lake (in m 2 ). This relationship has been used and proved reasonable for the approximation of glacial lake volume.

Streamflow
The study area is the origin of the Dudh Koshi River, which drains from the Everest region. The snow and glacier melt of the Everest region contribute a significant portion of the streamflow of the Dudh Koshi River [49]. The streamflow has been monitored at Rabuwa Bazar gauging station and made available by DHM, Nepal. Although WMO has recommended 1000 km 2 per station, this is the nearest available hydrological station representing the study area. In this study, the time-series data  were obtained from the records of the DHM. The outlier data were sliced by the IQR method. There were 5.56% of missing values in the monthly data. The monthly missing data were then replaced by the all-time (1995-2015) mean value of same month hydrological data. Thereafter, the homogeneity of data was examined using the RHTests software package. The monthly data were converted into four seasonal series (winter (DJF), spring (MAM), summer (JJAS), and autumn (ON)) and annual series by the averaging method. The seasonal and annual data from 2016 to 2019 were interpolated by a seasonal autoregressive integrated moving average (SARIMA) model. The model is written as SARIMA (p, d, q) × (P, D, Q, s), where p and seasonal P denote the number of time lags, d and seasonal D denote the degree of differencing, q and seasonal Q denote the number of moving average terms, and s denotes seasonal length in data [50]. The SARIMA (0, 1, 1) (0, 1, 1, 7) model was identified based on lowest Akaike information criterion (AIC) value, an estimation of prediction error. The trends of seasonal and annual time series were checked by the modified Mann-Kendall test and linear regression test.

Spatial Change Detection
The spatial change detection of alpine vegetation (1995 and 2020) was based on satellite image analysis. Landsat Surface Reflectance Datasets by USGS were used for change detection of spatial extent. These datasets were processed and atmospherically corrected data, which are available in the Google Earth Engine (GEE) dataset. The satellite images for the autumn period (September to December) were used for vegetation analysis. In general, in the annual cycle of vegetation growth, spring (MAM) is the start of the growing season for grass, shrub, and trees; summer (JJA) is the major growing season; autumn (SON) is the start of leaf-falling for deciduous plants; winter (DJF) is the dormant season for plants. However, alpine vegetations are mostly perennial. Hence, the autumn season is assumed to be suitable for vegetation monitoring by satellite image analysis. Since the output was only targeted at vegetation, image classification was done by the NDVI thresholding method. NDVI is the most widely used index in vegetation analysis. This was computed as the difference between near-infrared (NIR) and red reflectance values divided by their sum. NDVI value ranges from −1 to +1. NDVI values between 0.2 and 0.5 were classified as light vegetation, while NDVI values between 0.5 and 1 were classified as dense vegetation. Light vegetation in the Everest region includes alpine meadows/rangeland and juniper scrubland. Dense vegetation includes conifer forest and birch-rhododendron forest. To validate the classified outputs, 75 points for dense vegetation and 75 points for light vegetation were created in QGIS software ( Figure 3). Dense vegetation points were verified by ground observation, whereas light vegetation points were verified using Google Earth images. According to Congalton and Green (2019), the minimum number of sample points for accuracy assessment should be at least 20-100 points for each class [51].

NDVI Change Detection
USGS Landsat-5, -7, and -8 Tier 1 Surface Reflectance data were used for NDVI timeseries analysis. NDVI data during 1995-2013 were acquired from Landsat-5 data, but there were no data in 2002 because of sensor error [52] and 2012 due to near failure of Landsat-5 [53]. Although Landsat-7 has data from 1999 to the present, the failure of Scan Line Corrector of ETM+ sensor results in a zigzag pattern of image and loss of 22-25% of data, and hence they have limited use in scientific research [53]. The NDVI data from 2013 to 2019 were therefore extracted from Landsat-8.
The image collection was filtered according to regions of interest (ROI), i.e., alpine zone and period (1995-2020). Regarding quality assurance, these products contain pixel quality attribute information. The cloud and cloud shadow pixels were masked out. The NDVI of each image was calculated, and the mean NDVI (>0.2) of ROI was computed and extracted in Google Earth Engine (GEE). The unavailable NDVI data for 2002 and 2012 were filled by the mean value of corresponding monthly NDVI data. Further, the monthly NDVI data were converted into seasonal and annual series by the averaging method. The trend analysis was also done by using Mann-Kendall test and linear regression test.

Temperature Trend
The results of the modified Mann-Kendall test for annual and seasonal temperature data for 1995-2019 are presented in Table 2. A positive z-value (Table 2) and slope value ( Figure 4 and Table 3) indicate an increasing trend, whereas a negative z-value and slope value indicate a decreasing trend. The Kendall's tau value explains the strength of correlation between variables, ranging from −1 to +1. Annual and seasonal temperature trends based on the linear regression method are shown in Figure 4 and Table 3. Significances are provided at 90%, 95%, and 99% confidence levels.   In Table 2, the maximum annual temperature and minimum annual temperature trends show the existence of a monotonic trend at a 95% confidence level. Similarly, the mean annual temperature also shows the existence of a monotonic trend at a 90% confidence level. Besides, all the positive tau values and z-values in Table 2 indicate that the temperature is positively correlated or increasing with time series. As shown in Figure 4, variations from years to years are common in temperature trends but all indicate the warming trends. The maximum and minimum annual temperature trends are increasing significantly at a 95% confidence level, whereas the mean annual temperature trend shows a significant positive trend at a 90% confidence level. Similarly, the significant increasing trend is shown in the summer season of mean, maximum, and minimum temperature at a 99% confidence level (Table 2). According to the linear regression test shown in Figure 4, the mean annual temperature, maximum annual temperature, and minimum annual temperature show significant increasing trends. The magnitude of trends for mean, maximum, and minimum annual temperatures are 0.0329, 0.0415, and 0.055 • C/year, respectively.
In addition, summer maximum temperature and summer minimum temperature show significant increasing trends of 0.0343 and 0.0341 • C/year, respectively. The winter maximum temperature and winter minimum temperature show higher significant increasing trends of 0.078 and 0.0915 • C/year, respectively. This shows that winter warming is around 2.5 times faster than summer warming. Similarly, an insignificant but coherent warming trend is also found in spring. Autumn maximum temperature and autumn minimum temperature show warming trends of 0.0358 and 0.0467 • C/year, respectively.
Since ERA5 reanalysis data are interpolated data, they were again compared with the nearest observation data by the independent t-test. The results are shown in Table 4. We found that there is no significant difference between the mean of reanalysis and observation data.

Precipitation Trend
As indicated in Table 5, the annual precipitation shows a decreasing monotonic trend at a 95% confidence level. Winter and autumn precipitation show the existence of a monotonic trend at a 95% confidence level. With a 90% confidence level, summer precipitation also shows a negative monotonic trend. The negative z-value and tau value indicate the decreasing trend. In Figure 5, a significantly decreasing trend with a 90% confidence level is observed in the annual precipitation series. The magnitude of the trend is −14.90 mm/year.

Glacial Lake Changes
The results of the modified Mann-Kendall test for annual and seasonal data are shown in Table 5. It is observed that the annual surface water extent in the Everest region shows a positive monotonic trend at a 95% confidence level.
The seasonal and annual surface water extent trends based on the linear regression test are shown in Table 3 and Figure 6a. We found that the surface area of glacial lakes is increasing at the rate of 0.068 km 2 /year. The significant increasing trend is shown at a 99% confidence level. Based on the season, winter and autumn seasons show an increasing rate of 0.078 and 0.15 km 2 /year, respectively ( Table 3). The average surface area of glacier lakes in the study area is 4.9 km 2 . Thus, this observed increasing rate is quite significant. In Figure 6b, the glacial lake volume also shows a significant increasing trend with a rate of 0.0198 km 3 /year. The average glacial lake volume is 1.01 km 3 .

Streamflow Trend
According to the streamflow data for 1995-2015 recorded at Rabuwa Bazar gauging station, the maximum flowing rate is 2580 m 3 /s (cumecs), the minimum flowing rate is 9.07 m 3 /s, and the average flowing rate is 199 m 3 /s. The results of the modified Mann-Kendall test for annual and seasonal data are shown in Table 6 and Figure 7. In Table 6, it is observed that a significant decreasing trend is found in average annual streamflow and minimum annual streamflow at a 95% confidence level. The magnitudes of trend in average annual streamflow and minimum annual streamflow are −2.73 and −1.71 m 3 /s/year, respectively. It is also seen that the average summer streamflow and average minimum streamflow show significant decreasing trends. However, most of the trends for the maximum seasonal streamflow were not statistically significant. Similarly, all the spring season streamflows show significant decreasing trends.

Spatial Change Detection
Alpine vegetation represents life at climate limits and is a quick indicator of climate change. Based on GIS analysis, around 97.5% of the study area is a high-altitude region having an elevation of more than 3000 m. Thus, alpine vegetation in the Everest region mostly possesses alpine grasses, juniper shrubs, and conifer forests. Figure 8 shows the vegetation maps in 1995 and 2020. It is observed that the dense vegetation increased from 103.5 to 150 km 2 (45% increase), while light vegetation shows minor changes, increasing from 300.2 to 303.4 km 2 . This result shows the most marked expansion in dense vegetation in the Everest region. On the other hand, despite the constant areal extent, light vegetation such as alpine meadows/rangeland and juniper scrubs also shows expansion towards upper altitude, as shown in Figure 8. To check the classification accuracy, 75 random sample points for dense vegetation were examined by field observation and 75 random sample points for light vegetation were examined by Google Earth imagery (Figure 3). Then, the confusion matrix was created, which showed that the overall classification accuracy is 76.33% for the classified image. The classification accuracy for the dense vegetation is 89%, whereas the classification accuracy of light vegetation is 72% ( Table 7). The overall Kappa coefficient for 2020 is 0.76.

NDVI Change Detection
The result of the modified Mann-Kendall test for NDVI is shown in Table 5. It is revealed that there is a significant increasing trend in annual NDVI at a 95% confidence level. Similarly, all seasonal NDVI trends are also increasing significantly at a 95% confidence level. Figure 9 and Table 3 show the linear regression trend of annual and seasonal NDVI during 1995-2019. It is noted that the NDVI is increasing at a rate of 0.0035 units/year. Seasonally, the autumn NDVI shows the highest increasing rate of 0.0051 units/year. More importantly, the spring growing season also shows a significant increasing trend. An abrupt rise in annual NDVI was observed in 2008. This also indicates the greening trend is increasing in the Everest region.

Alpine Region Change and Its Relationship with Climate Variables
The responses of alpine vegetation (NDVI) and surface water (glacial lake area plus streamflow) are analyzed by the multiple linear regression method (Table 8). Here, the temperature and precipitation variables are considered as independent (x) variables, while the NDVI, surface water, and streamflow are considered as dependent (y) variables.
The table shows that 65% of the annual variability of NDVI is explained by mean, maximum, and minimum annual temperature. Further, regarding the seasonal variability, autumn temperature and summer precipitation explain 94% of the variability of annual NDVI. Similarly, 71% variability of surface water extent is represented by minimum annual temperature and annual precipitation. Overall, there is a high level of confidence that alpine vegetation and temperature are strongly correlated. Table 8. Multiple linear regression equations for the relationship between alpine vegetation (NDVI) and surface water (glacial lake plus streamflow) with temperature and precipitation.

Climate Trend
Temperature, precipitation, and solar radiation are the major climate controls that act as limiting factors for vegetation growth [54]. Some studies found that temperature is responsible for 31-33% of Earth's vegetation growth; precipitation is responsible for 40-52% of vegetation growth; solar radiation is responsible for 5-27% of vegetation growth [54,55]. In this study, it is revealed that the annual mean, maximum, and minimum temperatures at the Everest region are increasing at the rates of 0.0329, 0.0415, and 0.055 • C/year, respectively. In line with this result, Shrestha et al. [13] also found that the maximum temperature is increasing rate at 0.057 • C/year in the entire Nepal Himalayas. Within seasonal temperature, winter maximum temperature and winter minimum temperature show higher warming trends than the summer season. On the global scale, it is observed that winter warming is faster than summer warming [56]. Winter warming is ecologically more crucial. Winter warming can promote vegetation growth by reducing the seasonal frost of the winter season [57]. An experiment has discovered that soil respiration is increased by 9.3% by winter warming [58]. Winter warming stimulates litter decomposition in the winter season. Most importantly, prolonged winter warming triggers advancement in spring plant phenology and eventually increases the length of the growing season [59]. On the other hand, the annual precipitation of the Everest region shows a decreasing monotonic trend. However, the magnitude of change is found at −14.90 mm/year at a 90% confidence level. Besides this, summer precipitation also shows a decreasing trend at a 90% confidence level. Salerno et al. [60] also observed that the monsoon precipitation was weakening (−9.3 mm/year) in the Everest region during 1994-2013. Another larger-scale study including the Everest region by Shrestha et al. [61] found an increasing annual precipitation over the entire Koshi river basin (area of 87,970 km 2 ) from 1975 to 2010. However, this study did not include the stations in the study area, and the results in this study might be due to local spatial variation. Hence, the climate condition of the Everest region has changed over the recent 25 years period. This climate change is very likely to alter the alpine ecosystem and cryospheric processes of the Everest region.

Surface Water Dynamics
It is estimated that the glaciers are decreasing an average of 267 Gt/year at the global scale and 6 Gt/year in South Asia [62]. As a result, glacial lakes are growing. In the Everest region, the result shows the increasing trend of glacial lake surface area at 0.0676 km 2 /year. The glacial lake volume has been increasing at the rate of 0.0198 km 3 /year. Bajracharya et al. [63] also found that the moraine-dammed lakes' area, that is lakes formed due to blockade glacier debris, in the Dudh Koshi river basin increased from 2.291 to 7.254 km 2 between 1996 and 2000. Similarly, Salerno et al. [64] found that the number and area of glacier lakes in the Everest region especially in the designated area of Sagarmatha National Park have increased. They reported 624 glacial lakes with an area of 7.43 km 2 in 2008. The reason behind the glacial lake changes might be linked to temperature and precipitation change. A study of the Tibetan Plateau also found that the glacial lake expansion is highly correlated with glacier area recession, and this variation is associated with temperature rising [65]. The Everest region is also facing significant temperature rising, especially winter warming. Moreover, the glacial lake expansion may foster positive and negative impacts. There are unique glacier ecosystem services such as irrigation and drinking water supply to the downstream population and lowland discharge. However, most glacial lakes are formed by blockade of ice, glacier debris, landslide, and are inherently unstable type and susceptible to glacial lake outburst floods. Future endeavors should therefore be given to the study of glacial lake water resource utilization and their risk management.
Similarly, the average annual streamflow of the Dudh Koshi river basin is 199 m 3 /s. Dudh Koshi is a tributary of the Koshi River. The average flow of the Koshi River is 1658 m 3 /s. Similarly, the average flow of the three other major river basins of Nepal, namely, Gandaki, Karnali, and Mahakali, are 1753, 1441, and 698 m 3 /s, respectively [66]. In this study, from 1995 to 2019, we observed that the average annual streamflow is decreasing at the rate of 2.73 m 3 /s/year. These decreasing flow trends are consistent with weaker precipitation in the Everest region. Besides, the temperature is increasing, which might trigger more precipitation in the Everest region, but this is not the case. Gautam and Acharya [67] found a decreasing trend in the streamflow of the Dudh Koshi River basin from 1967 to 1997. Hence, the streamflow of the Dudh Koshi River is declining, which causes a water availability risk to the river-dependent ecosystems and societies. In addition to the need to closely monitor glacial lakes as mentioned above, it is also urgent and essential to evaluate the potential water stress associated with decreasing precipitation in this river basin.

Alpine Vegetation Dynamics
In this study, we found that NDVI is increasing at a rate of 0.00352 units annually. All the seasonal NDVIs are increasing with the highest in the autumn season. We also observed that annual NDVI is associated with autumn warming. Shi et al. found that autumn warning delays the end of the plant growing season for deciduous forests [68]. However, the study area consists of mostly evergreen conifer forests and perennial shrubs. Autumn warming also promotes vegetation growth by reducing frost damages [69]. The prominent increase in autumn NDVI reflects the increase in the extent of perennial vegetation. Additionally, the significant winter warming in the Everest region also favors vegetation growth. It is assessed that the dense vegetation cover has increased from 103.5 to 150 km 2 . Consistent with this, Baniya et al. [17] also observed a greening trend in all of Nepal with significant NDVI trends at 0.0005 and 0.002 units/year in the spring and autumn season, respectively. The driving factors for these increasing trends are attributed to temperature (39.62%), ecological restoration (25.16%), and multiple factors (35.22%). Over the larger scale in the Koshi river basin, Wu et al. [18] found a greening trend of 0.0023 units/year from 2000 to 2018. Zhang et al. [15] also showed the increasing rate of average growing season NDVI over the Koshi River basin at the rate of 0.0034 units/year from 2000 to 2011.
In the context of greening in the Everest region, there might be multiple casual/control factors. The greening trend is increasing despite weaker precipitation. The ecological restoration program was also held in the study area. However, the area of the plantation during 1979-2008 is small (1.51 km 2 ) as compared to the increasing area [9]. Ecological restoration in the Everest region is also challenging because of less water, cold weather, and rocky topography. Besides, the Everest region had around 1.02 km 2 settlement area and 8.50 km 2 agriculture area in 2011, and this region is managed under park regulations. Hence, the human disturbances to land-use change are considered minimal [28]. The shrub promotion due to overgrazing and abandonment of remote farmland is likely to favor vegetation expansion. Hence, it is obvious that the vegetation expansion is also driven by climatic factors. In this study, we found that annual NDVI is highly correlated with annual temperature. The abrupt rise of annual NDVI during 2008-2015 coincided with an increasing trendline in annual mean temperature during 2003-2009. Moreover, the annual precipitation showed an increasing trendline during 2008-2012 despite weak precipitation over the period. Furthermore, the temperature change can drive an upper shift in snow line altitude, which subsequently favors temperature for vegetation expansion and migration [70]. It is also revealed that temperature triggered alpine vegetation greening on the Tibetan Plateau from 1982 to 2011 [71]. For future endeavors, long-term research on phenological changes and carbon dynamics of the alpine vegetation may explore more about the climate change impacts on the high-altitude region.

Conclusions
This study demonstrates the changes in climate variables and their relationships with the change extents in surface water and alpine vegetation of the Everest region from 1995 to 2019. It is revealed that the mean, maximum, and minimum annual temperature increasing rates are 0.0329, 0.0415, and 0.055 • C/year, respectively. These increasing rates are statistically significant at a 95% confidence level. On the other hand, the annual precipitation also shows a significantly decreasing monotonic trend. The surface water area of glacial lakes shows an increasing rate of 0.0676 km 2 /year, whereas the glacial lake volume shows an increasing trend of 0.0198 km 3 /year. On the other hand, the average annual streamflow of the Dudh Koshi river basin is also decreasing at the rate of 2.73 m 3 /s/year. We also found that the alpine vegetation extent of the Everest region is changing. The dense vegetation cover in 2020 is 150 km 2 , a 46.5 km 2 increase from the vegetation area in 1995. The greening trend in the Everest region is significant as the NDVI is increasing at a rate of 0.0035 units annually. In the context of seasonal variability, winter warming is 2.5 times faster than summer warming. The glacial lake change in the winter is increasing at the rate of 0.078 km 2 /year. The most prominent NDVI increase in the autumn season indicates vegetation expansion. Regarding the correlation of climate variables and NDVI, surface water extent, and streamflow, we found that annual NDVI is significantly correlated with maximum annual temperature, mean annual temperature, and minimum annual temperature. We also found that glacial lake expansion shows a significant association with minimum annual temperature and annual precipitation. Overall, this study reveals a significant alteration in the alpine ecosystem of the Everest region that could impact on the water-energy-food nexus of the Dudh Koshi river basin. The results indicate an urgent need for evaluation of the climate change impacts and implementation of effective planning to cope with threats from glacial melting and decreasing trends in precipitation, along with the expansion of alpine vegetation.