Spatial and Temporal Variation of NDVI in Response to Climate Change and the Implication for Carbon Dynamics in Nepal

Nepal is a country of contrast, with varying altitude, climate and vegetation from the top of high mountains to the tropical forest in low lands. The terrestrial vegetation has rapidly been altered by climate change in Nepal. The spatial and temporal evolution of vegetation and its linkage to climatic variables were analyzed using the Normalized Difference Vegetation Index (NDVI) obtained from Advanced Very High Resolution Radiometer (AVHRR) sensors. A linear regression model and Sen’s slope method were used to estimate NDVI trends and the Pearson correlation between NDVI and climatic variable, i.e., temperature and precipitation were calculated to identify the role of climate in vegetation changes. The carbon dynamics were also measured using a biomass carbon density estimation model. Results showed that NDVI experienced an overall increasing trend in Nepal from 1982–2015. The NDVI significantly increased at the rate of 0.0008 year−1 (p < 0.05) with seasonal variation of 0.0004 year−1, p > 0.05; 0.0007 year−1, p < 0.05; 0.0008 year−1, p < 0.05 and 0.0007 year−1, p > 0.05 in winter, pre-monsoon, monsoon and post-monsoon seasons, respectively. The NDVI relative change ratio (RCR) was 6.29% during last 34 years in Nepal. The correlation between NDVI and temperature was significantly positive (r = 0.36, p = 0.03), but there was a negative correlation with precipitation (r = −0.21, p = 0.28). Altogether, 82.20% of the study areas showed a positive correlation with temperature in which 34.97% was significant and 69.23% of the area had a negative correlation (16.35% significant, p < 0.05) with precipitation. In addition, NDVI-based carbon estimation showed that Nepal’s forest total carbon stock is 685.45 × 106 t C (i.e., an average of 115.392 t C/ha) with an annual carbon sequestration rate of 0.10 t C/ha from 1982–2015. The results suggest that NDVI variation is more sensitive to temperature than precipitation and it is valuable to measure carbon dynamics in Nepal.


Introduction
Research on the vegetation response to a global change has been a central activity in earth system science during the past two decades [1].Vegetation is strongly influenced by climate change, i.e., temperature, precipitation and radiation [2][3][4], which are also driven by elevated atmospheric CO 2 concentration and changing rates of nitrogen deposition [5].Climate change has increased the potential risk in ecosystems in several parts of the earth [6] and has rapidly altered the earth system dynamics [7].The monitoring of these vegetation changes using earth observing satellite data has improved our understanding of vegetation dynamics [8,9].The satellite-derived Normalized Difference Vegetation Index (NDVI), which is the normalized ratio of red and near-infrared (NIR) reflectance [10], can be regarded as an indicator of vegetation activities [1].The NDVI has been used to detect vegetation dynamics at regional and global scales [11][12][13][14].It can also be used to determine biotic responses to climate change [2].A positive NDVI trend in global land ecosystem has been reported based on the remarkable 30-year archive of satellite imagery [15].Results of previous studies suggest increasing vegetation in the Northern Hemisphere [1,[16][17][18][19][20][21], which is generally controlled by temperature rise in middle and high latitudes [22].The NDVI can be used to measure ecological variables like plant productivity, photosynthetic activity, plant phenology, tropic interactions, biomass and the global carbon cycle [23].The relationship between NDVI and climate variables has also been extensively studied worldwide [14,22,24].These previous studies suggest that temperature and precipitation are two major climatic factors that are closely linked to vegetation activities by interrupting the photosynthetic process [21].The NDVI is positively correlated with temperature and negatively correlated with precipitation in the Northern Hemisphere [21].Vegetation plays a leading role in global carbon cycles [25] which can be measured using NDVI through an accumulation difference between gain (photosynthesis) and loss (respiration).The NDVI has a multidimensional implication and has also been used for land cover classification [26], drought monitoring [27,28], land degradation [29] and forest carbon dynamics [30][31][32].
In Nepal, vegetation is sensitive to climate change [33] because it is an important indicator of ecosystem dynamics [34].The abrupt changes in the climate, i.e., increased temperature and rainfall variability in high topographic regions of Nepal and surrounding Himalaya, resulted in an upward shifting of vegetation, the presence of new invasive species, prolonged dry spells, and higher incidence of disease and forest fire, which can alter structure and function, productivity and seasonal phenology of the vegetation.The temperature showed an increasing trend at the rate of 0.06 • C year −1 in most of the hilly areas and mountains from 1971-1994 [35].Precipitation presents temporal and spatial variations in Nepal [36].The summer monsoon dominates the total precipitation whereas the winter monsoon contributes only 3% of total annual precipitation [37].Therefore, temperature and monsoon variation has a profound impact on vegetation growth.Several studies suggest that Nepal is one of the more vulnerable countries with a rapidly altering ecosystem and land surface processes due to the influence of climate change [38].Heavy rainfall events, natural disasters and severity of drought are causing problems in the structure and function of the ecosystem [39].The forest, biodiversity and natural ecosystems are main thematic sectors which are more severely impacted by climate change [33].The large population growth [40], cropland expansion [41], deforestation and land uses changes have unevenly interrupted the vegetation system in different geographical regions [42].This tends to suggests that rapid alteration of the vegetation can cause widespread damage in terrestrial ecosystems, agriculture, bio-diversity and water resources [43].There are big challenges to address these bio-physical environmental problems by managing good ecological systems.Therefore, understanding vegetation dynamics is crucial for better environment management in Nepal.There are numerous studies on vegetation dynamics in this region, but the use of a remote sensing approach is limited.Recently, one of the studies showed positive NDVI trends which were mainly attributed to CO 2 fertilization in Nepal [44].In the tree line ecotone, most of the studies were related to tree rings.The relationship between tree rings and climate change has been extensively studied in Nepal [45][46][47][48][49][50][51][52] but these studies were site specific, species specific and used a single approach of dendrochronology.Most of the research followed a similar pattern by examining tree ring growth, identifying regeneration and analyzing quantitative characteristics.However, the relationship between vegetation and climatic changes is poorly understood at large temporal and spatial scales.A ground-based inventory was conducted to estimate forest carbon stocks in Nepal [53] which only considered those trees with a diameter greater than 10 cm at breast height.In this context, a satellite-derived vegetation index could be used as a surrogate to study vegetation changes and carbon dynamics in Nepal.
This study aims to gain a better understanding of the changing spatial and temporal pattern of NDVI in response to varying temperature and precipitation and its implication on carbon dynamics in Nepal over the period of 1982-2015.Specifically, the study was planned to identify long-term NDVI trends in Nepal so that the existing and changing paradigm of physical and biological environmental conditions can be determined for better ecology and environmental management.Second, the correlation between NDVI and climatic factors was identified to determine the role of climate in vegetation change at large temporal and spatial scales.Identification of long-term climate and vegetation change itself is crucial to deal with the alteration of livelihood, ecosystem dynamics and land-atmosphere interactions.Knowledge of vegetation and climate changes in different ecological regions has significant influences on ecological systems in Nepal because the Terai, Hills and Mountain areas represent different altitudes, ecological patterns and contrasting bio-physical environment [54].Finally, the NDVI was used to identify the total forest carbon stock and carbon trends in Nepal; this is an indicator of ecosystem dynamics, indicating whether the ecosystem is restoring or degrading.For our purpose, we used the most popular NDVI data [55] due to their long and continuous record [56].These have been widely used because of the high quality and reduced noise from volcanic eruptions, solar angles and sensor errors [2,57,58].These data are considered as the longest temporal coverage and most suitable data for analysis of vegetation dynamics [59].This data set has also proved to be one of the best products to detect temporal changes in vegetation [60].In this study, we examined annual and seasonal NDVI changes, NDVI trends in low, middle and high lands of different altitudinal zones and the relationship of NDVI changes with temperature and precipitation.In addition, the spatial and temporal variation of forest carbon was measured based on the NDVI approach.Remote sensing-based carbon estimation is a new approach in Nepal which can be practiced at large temporal and spatial scales.This research was designed at a landscape scale for the entire country instead of point-, localand site-specific purposes, which are generally used for traditional data at small temporal and spatial scales.In addition, the study aims to determine the relevancy of remote sensing in vegetation and climate change research.Here, we grouped the vegetation into forest, shrub, grass and crop lands [61].

Study Area
Nepal is located between north latitude 26 • 22 and 30 • 27 and east longitude 80 • 04 and 88 • 12 in south Asia.It is bounded by the Tibetan highland to the north and an Indian foothill of the Himalaya to the south.The altitude varies between 60-220 m in the south and reaches a maximum of 8848 m in the north.It covers an area of 147,181 km 2 with a length of about 800 km, parallel to the Himalayan axis, and a width of 160 km.There are complex ecosystems, and the bio-climatic zone extends from the broad leaved forest in the tropical climatic zone below 1000 m altitude to an alpine forest characterized by various stunted bushy shrubs above 4000 m altitude.The majority of the land is forest, which occupied 39.1% of the total area [62].This value increased to 40.36% of the total area during 2015 [53].
The climatic gradient of Nepal varies due to complex topography and different physiographic features.There are three main ecological zones, i.e., Terai (low lands), Hills (midlands) and Mountains (high lands) in Nepal (Figure 1): the elevation ranges below 1000 m altitude in Terai, from 1000-3000 m altitude in the hills and an altitude greater than 3000 m represents high lands, i.e., the mountains and Himalaya [54].The Terai, Hill and Mountain areas occupy 21, 36 and 20 districts, respectively, in Nepal [63].

NOAA, NDVI3g.v1 and Climate Data Sets
The NDVI data set used in this study is the latest version of the dataset, namely, NDVI3g.v1downloaded from https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1 [64].The data were generated by NOAA's Advanced Very High Resolution Radiometer (AVHRR) sensors under the framework of the Global Inventory Monitoring and Modeling System (GIMMS).The data were processed by the GIMMS team using adaptive empirical model decomposition [65] to address the issues introduced by calibration, sensor degradation, view geometry, orbital drift, atmospheric contamination and sampling related to the AVHRR sensor.The GIMMS NDVI product has an 8 km × 8 km spatial resolution and a 15-day temporal resolution.The data set covers the period from 1982-2015.In this study, the temperature and precipitation datasets were obtained from the latest version of the Climate Research Unit (CRU) time series 4.1 data sets from the University of East Anglia [66][67][68].These are gridded data sets with a spatial resolution of 0.5° for the period 1982-2015.

Methods
The NDVI and its change attributed to climate change were investigated based on NDVI, temperature and precipitation time series data.Altogether, 814 NDVI images (24 image/year) for the periods of 34 years were downloaded, and raster images were analyzed.The monthly NDVI data were generated from bi-monthly data using the maximum value composite method [69].In this study, the pixels with annual mean NDVI values less than 0.1 were removed [2].The annual, seasonal and pixel-wise NDVI trends were computed.The seasons were mainly categorized as winter (December-February), pre-monsoon (March-May), Monsoon (June-September) and post-monsoon (October-November) based on temperature and precipitation patterns in Nepal.We downscaled CRU data from 0.5 degrees to the spatial resolution of 8 km to ensure similar resolution with NDVI data using Delta and Hermite Interpolating Polynomial (PCHIP) interpolation methods [70].The NDVI trends in different ecological zones were calculated using the data layer of the ecological zone obtained from the Department of Survey and raster NDVI of the mean, RCR and slope using raster calculation tools in ArcGIS.

NOAA, NDVI3g.v1 and Climate Data Sets
The NDVI data set used in this study is the latest version of the dataset, namely, NDVI3g.v1downloaded from https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1 [64].The data were generated by NOAA's Advanced Very High Resolution Radiometer (AVHRR) sensors under the framework of the Global Inventory Monitoring and Modeling System (GIMMS).The data were processed by the GIMMS team using adaptive empirical model decomposition [65] to address the issues introduced by calibration, sensor degradation, view geometry, orbital drift, atmospheric contamination and sampling related to the AVHRR sensor.The GIMMS NDVI product has an 8 km × 8 km spatial resolution and a 15-day temporal resolution.The data set covers the period from 1982-2015.In this study, the temperature and precipitation datasets were obtained from the latest version of the Climate Research Unit (CRU) time series 4.1 data sets from the University of East Anglia [66][67][68].These are gridded data sets with a spatial resolution of 0.5 • for the period 1982-2015.

Methods
The NDVI and its change attributed to climate change were investigated based on NDVI, temperature and precipitation time series data.Altogether, 814 NDVI images (24 image/year) for the periods of 34 years were downloaded, and raster images were analyzed.The monthly NDVI data were generated from bi-monthly data using the maximum value composite method [69].In this study, the pixels with annual mean NDVI values less than 0.1 were removed [2].The annual, seasonal and pixel-wise NDVI trends were computed.The seasons were mainly categorized as winter (December-February), pre-monsoon (March-May), Monsoon (June-September) and post-monsoon (October-November) based on temperature and precipitation patterns in Nepal.We downscaled CRU data from 0.5 degrees to the spatial resolution of 8 km to ensure similar resolution with NDVI data using Delta and Hermite Interpolating Polynomial (PCHIP) interpolation methods [70].The NDVI trends in different ecological zones were calculated using the data layer of the ecological zone obtained from the Department of Survey and raster NDVI of the mean, RCR and slope using raster calculation tools in ArcGIS.

Linear Regression Method (LRM)
The linear regression method was applied to analyze annual and seasonal NDVI trends [71].The slope estimated by LRM indicated the mean temporal change of NDVI (Equation ( 1)).A positive slope indicates increasing trends while negative slopes indicate decreasing trends [72,73].
where slope is the trend of the vegetation dynamics, n is the number of years in the study period, i is the year and NDVI i is the mean NDVI in the ith year.A t-test [74] in MATLAB was used for each pixel for significant trend analysis of the NDVI slope.

Sen's Non-Parametric Slope
In addition to the linear regression method, Sen's slope was also used to estimate the slope of NDVI in Nepal.If the time series data present a linear trend, the true slope (change per unit time) of a trend could be estimated by the non-parametric index developed by Sen [75], which is based on the assumption of a linear trend.In the last three decades, it has been highly applied in the studies of vegetation dynamics [14,[76][77][78]].Sen's slope was computed as where x i and x j are the NDVI values at time i and j, respectively (i > j).The median is the median function.The negative values of Sen's slope indicate a negative trend and a positive value indicates a positive trend.The significance of vegetation changes was evaluated using the nonparametric Mann-Kendall test [79,80] with statistical significance at the confidence interval of 95% (p < 0.05).
H o (no trend) was rejected if the p value was less than a predefined significance level of 0.05.The Kendall tau (−1 to +1) obtained from Mann-Kendell statistics was also used to measure the strength of the relationship between two variables.

Relative Change Ratio (RCR)
The annual and seasonal relative change ratios (RCR) that either increase or decrease were estimated using the slope and mean NDVI values [77].The RCR of NDVI, which is also called greenery change ratio (%), was used during the study period.
where slope is the slope of annual NDVI estimated by linear regression or Sen's slope (here we used the linear regression slope), mean is the average annual NDVI values and N is the length of the study period from 1982-2015.

Correlation between NDVI and Climatic Variables
The Pearson correlation coefficient (r) between NDVI and climatic factors was calculated to assess the impact of climate change on vegetation dynamics.This approach has been widely applied to analyze the relationship between NDVI and climatic factors [14,[81][82][83].The r indicates the linear relationship and the p values indicate the significant correlation level between two variables.The Pearson correlation for each pixel was calculated using a statistical package in R, and a t-test [74] was performed for significant trend analysis.If the correlation value between two variables is positive and p value is less than 0.05, the correlation is statistically significant.Here, the correlation map was developed using annual raster data sets of NDVI, temperature and precipitation.

Biomass Carbon Density (BCD) Estimation Model
This model is the satellite-based empirical global forest biomass model used by [32] to estimate biomass carbon density (BCD) based on corresponding values of NDVI max , latitude and longitude (Equation ( 4)).This model was also used by [30,31] based on growing season total NDVI, annual average NDVI and latitude.BCD = 111.521ln(NDVI) − 0.452 lat − 20.034 lon + 0.08568 lon 2 + 1278.29 (4) where BCD is forest biomass C density (Mg C/ha), NDVI is the NDVI max , and lat and long are the latitude and longitude.In our approach, we also used a corresponding value of NDVI max as one of the predictors of biomass carbon density because it is difficult to clearly define the length of the growing season in a coniferous forest, and we assumed that NDVI max values correspond to the optimal measurement condition [69].Therefore, the input variables in our study were NDVI max , latitude and longitude of the area.Here, the forest biomass refers to the total above and below ground biomass.

Temporal NDVI Trends
The NDVI trends were positive in estimation of both the linear and Sen's slope.Based on Sen's slope, the annual NDVI has significantly increased at the rate of 0.0008 year −1 .The mean and standard deviation of the annual NDVI from 1982-2015 was 0.50 and ±0.015, respectively.The relative change ratio of NDVI was 6.29% during the study period in Nepal (Table 1).The NDVI trend was positive in all the seasons.The pre-monsoon and monsoon season NDVI has significantly increased at the rate of 0.0007 year −1 (R 2 = 0.10, p = 0.03) and 0.0008 year −1 (R 2 = 0.23, p = 0.01), respectively (Table 1; Figure 2b,c) while NDVI showed positive but insignificant trends during winter and post-monsoon seasons (Table 1; Figure 3b,e).In winter, the mean and standard deviation of the NDVI were 0.491 and ±0.015, respectively.The NDVI increased at the rate of 0.0006 year −1 in the winter season (Sen's slope-0.0004year −1 , p = 0.1).The RCR in winter is 4.14%, which is lower than that of other seasons, i.e., 4.89%, 5.93% and 5.41% in pre-monsoon, monsoon and post-monsoon seasons, respectively.Both the linear and Sen's slope for annual and seasonal time scales are given in Table 1.The mean NDVI in winter was relatively higher than the pre-monsoon and lower than the annual, monsoon and post-monsoon seasons.In the pre-monsoon season, the mean NDVI was 0.43, which was lower than that of the annual, monsoon and post-monsoon seasons.The mean NDVI values were 0.50 and 0.58 in the monsoon and post-monsoon seasons, respectively; these were the highest mean NDVI values compared to other seasons.The NDVI increased with a trend of 0.0007 year −1 in pre-monsoon and post-monsoon seasons.The annual and seasonal standard deviations were low (Table 1), and the annual Kendall tau of NDVI was 0.45 from 1982-2015.The positive Kendall's tau statistics also indicate a positive NDVI trend.The figures above were developed based on obtained Mann-Kendall test statistics and Sen's slope values which showed annual and seasonal NDVI trends and inter-annual variation of mean NDVI distribution in Nepal from 1982 to 2015.The analysis of temporal NDVI changes in Terai (<1000 m), Hill (1000-3000 m) and Mountain (>3000 m) areas showed that NDVI has highest in the hills, with an annual positive trend of 0.0012, followed by Terai with a mean NDVI of 0.56 and annual average trend of 0.001 (Table 2).The mountain area, i.e., above 3000 m altitude, has a low mean NDVI and annual trends.In general, the hilly parts of Nepal have high mean NDVI and increasing trends at annual and seasonal time scales compared to the mountains and the Terai (Table 2).The mean NDVI was higher during The mean NDVI in winter was relatively higher than the pre-monsoon and lower than the annual, monsoon and post-monsoon seasons.In the pre-monsoon season, the mean NDVI was 0.43, which was lower than that of the annual, monsoon and post-monsoon seasons.The mean NDVI values were 0.50 and 0.58 in the monsoon and post-monsoon seasons, respectively; these were the highest mean NDVI values compared to other seasons.The NDVI increased with a trend of 0.0007 year −1 in pre-monsoon and post-monsoon seasons.The annual and seasonal standard deviations were low (Table 1) 2).The mountain area, i.e., above 3000 m altitude, has a low mean NDVI and annual trends.In general, the hilly parts of Nepal have high mean NDVI and increasing trends at annual and seasonal time scales compared to the mountains and the Terai (Table 2).The mean NDVI was higher during the post-monsoon season in all ecological zones in Nepal.The NDVI RCR and trends were also positive in all ecological zones, with smaller changes shown for the mountain.The NDVI RCR over the period of 34 years was 7.20% in the hills, followed by 7.16% in Terai.This showed that mountains have a low mean NDVI value but positive NDVI trends in Nepal.

Spatial Variation of NDVI
The mean NDVI and trends were spatially varied in different ecological regions.The annual mean NDVI in northern parts was very low with negative NDVI in some pixels, and the NDVI value was high (above 0.4) in the middle zone of Nepal (Figure 3).The NDVI generally ranged from 0.  the post-monsoon season in all ecological zones in Nepal.The NDVI RCR and trends were also positive in all ecological zones, with smaller changes shown for the mountain.The NDVI RCR over the period of 34 years was 7.20% in the hills, followed by 7.16% in Terai.This showed that mountains have a low mean NDVI value but positive NDVI trends in Nepal.

Spatial Variation of NDVI
The mean NDVI and trends were spatially varied in different ecological regions.The annual mean NDVI in northern parts was very low with negative NDVI in some pixels, and the NDVI value was high (above 0.4) in the middle zone of Nepal (Figure 3).The NDVI generally ranged from 0.4 to 0.5 in the low lands of southern regions.The maximum mean NDVI ranges from 0.7 to 0.8, mainly in mid-hill areas occupied by broadleaved dense forest.The average NDVI values in western Nepal (from the western boundary to 83° E) and central Nepal (from 83° to 86°3′ E) were 0.481 and 0.485, respectively, but 0.551 in eastern Nepal (from 86°3′ E to eastern boundary).The NDVI in 83.89% (57.35% significant, p < 0.05) of the study area showed an increasing trend, and 16.10% (4.68% significant, p < 0.05) of the area presents a decreasing trend (Figure 4a).The red colors in the small figures represent significant changes, where the p value is less than 0.05, and the green colors represent insignificant changes (p > 0.05).The spatial variations of the NDVI trends were seasonally different, with a range from −0.01 to 0.01 year −1 .The annual NDVI trends were negative, mainly in the north-west and some parts of southern regions, with a maximum range of 0.01 year −1 .Similar patterns were identified in the north-west region during the post-monsoon season (Figure 4e).In the pre-monsoon season, southern and eastern parts experienced negative trends in some pixels.In winter, 72.73% of the area showed a positive trend, which is higher than in the pre-monsoon season (68.41%) but lower than in the post-monsoon season (88.63%).The NDVI in 83.89% (57.35% significant, p < 0.05) of the study area showed an increasing trend, and 16.10% (4.68% significant, p < 0.05) of the area presents a decreasing trend (Figure 4a).The red colors in the small figures represent significant changes, where the p value is less than 0.05, and the green colors represent insignificant changes (p > 0.05).The spatial variations of the NDVI trends were seasonally different, with a range from −0.01 to 0.01 year −1 .The annual NDVI trends were negative, mainly in the north-west and some parts of southern regions, with a maximum range of 0.01 year −1 .Similar patterns were identified in the north-west region during the post-monsoon season (Figure 4e).In the pre-monsoon season, southern and eastern parts experienced negative trends in some pixels.In winter, 72.73% of the area showed a positive trend, which is higher than in the pre-monsoon season (68.41%) but lower than in the post-monsoon season (88.63%).A high proportion of the land has a positive trend during the post-monsoon season, whereas only 11.36% of the area showed a negative NDVI trend, especially the north-western parts of Nepal.The NDVI has significantly changed in 46.39%, 38.88%, 43.46% and 40.99% of the study area in the winter, pre-monsoon, monsoon and post-monsoon seasons, respectively.Most had positive significant trends.An annual significant trend was mostly found in middle hills and south-eastern parts of Nepal where the land is mostly occupied by broad leaved and needle leaved open and closed forests and agricultural land (Figure 1).The NDVI change has seasonally and spatially varied in which western parts had more significant NDVI trends during both winter and monsoon seasons.Eastern Terai has significant NDVI trends at annual and all seasonal time scales.Some of the pixels also had significant negative trends during the pre-monsoon season compared to others.The large fraction of the significant trends showed high vegetation dynamics during the last three decades in Nepal.

Relationship between NDVI and Climate Change
The annual NDVI and temperature have significantly increased but the precipitation decreased (NDVI: 0.0009 year −1 , temperature: 0.03 °C year −1 and precipitation: −5.12 mm year −1 ).The annual trend of the precipitation was negative, but the annual NDVI trend was positive.The temperature in pre-monsoon and post-monsoon seasons increased by 0.04 °C year −1 , precipitation decreased by 0.47 mm year −1 and 3.94 mm year −1 but the NDVI increased by 0.0006 year −1 and 0.0007 year −1 , respectively.The precipitation trend was more negative during the monsoon season (3.94 mm year −1 ) where the temperature increased at the rate of 0.03 °C, and NDVI also increased significantly.The Pearson correlation analysis showed that the correlation between annual NDVI and temperature was significantly positive (r = 0.36, p = 0.03) except during the pre-monsoon season (r = −0.05,R 2 = 0.005).A high proportion of the land has a positive trend during the post-monsoon season, whereas only 11.36% of the area showed a negative NDVI trend, especially the north-western parts of Nepal.The NDVI has significantly changed in 46.39%, 38.88%, 43.46% and 40.99% of the study area in the winter, pre-monsoon, monsoon and post-monsoon seasons, respectively.Most had positive significant trends.An annual significant trend was mostly found in middle hills and south-eastern parts of Nepal where the land is mostly occupied by broad leaved and needle leaved open and closed forests and agricultural land (Figure 1).The NDVI change has seasonally and spatially varied in which western parts had more significant NDVI trends during both winter and monsoon seasons.Eastern Terai has significant NDVI trends at annual and all seasonal time scales.Some of the pixels also had significant negative trends during the pre-monsoon season compared to others.The large fraction of the significant trends showed high vegetation dynamics during the last three decades in Nepal.

Relationship between NDVI and Climate Change
The annual NDVI and temperature have significantly increased but the precipitation decreased (NDVI: 0.0009 year −1 , temperature: 0.03 • C year −1 and precipitation: −5.12 mm year −1 ).The annual trend of the precipitation was negative, but the annual NDVI trend was positive.The temperature in pre-monsoon and post-monsoon seasons increased by 0.04 • C year −1 , precipitation decreased by 0.47 mm year −1 and 3.94 mm year −1 but the NDVI increased by 0.0006 year −1 and 0.0007 year −1 , respectively.The precipitation trend was more negative during the monsoon season (3.94 mm year −1 ) where the temperature increased at the rate of 0.03 • C, and NDVI also increased significantly.The Pearson correlation analysis showed that the correlation between annual NDVI and temperature was significantly positive (r = 0.36, p = 0.03) except during the pre-monsoon season (r = −0.05,R 2 = 0.005).
The seasonal correlation between NDVI and precipitation was negative except in winter and post-monsoon seasons where the correlation was slightly positive.In the pre-monsoon season, the correlation between NDVI and temperature was negative.The NDVI has a positive correlation with both temperature and precipitation in winter and post-monsoon seasons (Table 3).The pixel-wise correlation between NDVI and annual average temperature and precipitation from 1982-2015 is shown in Figure 5.The seasonal correlation between NDVI and precipitation was negative except in winter and post-monsoon seasons where the correlation was slightly positive.In the pre-monsoon season, the correlation between NDVI and temperature was negative.The NDVI has a positive correlation with both temperature and precipitation in winter and post-monsoon seasons (Table 3).The pixel-wise correlation between NDVI and annual average temperature and precipitation from 1982-2015 is shown in Figure 5.There were positive correlations between annual average temperature and NDVI in 82.20% of the study area (34.97% is significantly positive, p < 0.05).However, there were negative correlations in 17.79% of the study area in which 2.10% was significant (p < 0.05).Similarly, 30.76% of the study area has a positive correlation between NDVI and precipitation, with 1.08% significantly positive (p < 0.05), but 69.23% of the area has a negative correlation in which 16.35% is significantly negative (p > 0.05).The central hills have high positive and significant correlations in the majority of the area between NDVI and temperature, with a maximum correlation coefficient of 0.73.The negative correlation was mostly observed in eastern, southern and northern Nepal with a maximum negative correlation of 0.59 between NDVI and temperature (Figure 5a).High negative correlations between NDVI and precipitation were mainly observed in south-western region.The majority of this area showed significant negative correlations with between NDVI and precipitation (Figure 5d).The maximum positive correlation was 0.46 and the negative correlation was 0.60 between NDVI and precipitation.The significant positive correlation was higher between NDVI and temperature compared to the NDVI and precipitation.There were very few areas with significant positive There were positive correlations between annual average temperature and NDVI in 82.20% of the study area (34.97% is significantly positive, p < 0.05).However, there were negative correlations in 17.79% of the study area in which 2.10% was significant (p < 0.05).Similarly, 30.76% of the study area has a positive correlation between NDVI and precipitation, with 1.08% significantly positive (p < 0.05), but 69.23% of the area has a negative correlation in which 16.35% is significantly negative (p > 0.05).The central hills have high positive and significant correlations in the majority of the area between NDVI and temperature, with a maximum correlation coefficient of 0.73.The negative correlation was mostly observed in eastern, southern and northern Nepal with a maximum negative correlation of 0.59 between NDVI and temperature (Figure 5a).High negative correlations between NDVI and precipitation were mainly observed in south-western region.The majority of this area showed significant negative correlations with between NDVI and precipitation (Figure 5d).The maximum positive correlation was 0.46 and the negative correlation was 0.60 between NDVI and precipitation.The significant positive correlation was higher between NDVI and temperature compared to the NDVI and precipitation.There were very few areas with significant positive correlation between NDVI and precipitation.This pixel-wise correlation showed that NDVI generally has a positive correlation with temperature and a negative correlation with precipitation.

Implication of NDVI for Carbon Dynamics in Nepal
The average forest carbon stock in Nepal was 115.392 t C/ha with an annual carbon sequestration rate of 0.10 t C/ha from 1982-2015.Nepal's total forest carbon stock was 685.45 × 10 6 t C. The average annual carbon stock was high in the years of 1990, 1995 and 2007 and low in 1982 and 1985 (Figure 6a).It can be concluded that the carbon stock was lower before the 1990s and it gradually increased thereafter.
Forests 2018, 9, x FOR PEER REVIEW 11 of 18 correlation between NDVI and precipitation.This pixel-wise correlation showed that NDVI generally has a positive correlation with temperature and a negative correlation with precipitation.

Implication of NDVI for Carbon Dynamics in Nepal
The average forest carbon stock in Nepal was 115.392 t C/ha with an annual carbon sequestration rate of 0.10 t C/ha from 1982-2015.Nepal's total forest carbon stock was 685.45 × 10 6 t C. The average annual carbon stock was high in the years of 1990, 1995 and 2007 and low in 1982 and 1985 (Figure 6a).It can be concluded that the carbon stock was lower before the 1990s and it gradually increased thereafter.The forest biomass carbon density showed spatial heterogeneity (Figure 6b).It was high in the south-western and central areas but lower in eastern and northern parts of Nepal.The highest carbon density reached 168.82 t C/ha in western parts and the lowest was 9.27 t C/ha in upper parts of hills mainly occupied by sub-alpine forests.In Figure 6b, the green color represents the dense broad leaved forest.The light green colors also represent the broad leaved open and closed forest with few needle leaved forest.Similarly, the red and yellow colors mostly represent needle leaved forests in Nepal (Figure 1).

Discussion
The NDVI experienced increasing trends in large areas of Nepal from 1982-2015 (Table 1; Figures 2 and 4).The annual NDVI increased by 0.0008 year −1 with positive trends in the estimation of both the linear and Sen's slope for all seasons, with obvious spatial heterogeneity during the last three decades.The majority of the area has positive NDVI trends with negative trends in a few locations.The results are consistent with previous studies carried out in several parts of the world, mainly in the Northern Hemisphere [17,19,21,78,84].Vegetation change is regarded as an indicator of environmental change which has severe ecological consequences in Nepal.The ecology of vegetation changes has significance impacts on functioning ecosystems, energy balance, species richness and land atmospheric interactions.Generally, high NDVI values (indicated by an integrated NDVI or the annual maximum NDVI value) reflect a good status of plant productivity and biomass.A positive NDVI trend indicates that ecology and ecosystem conditions are better.The ecological restoration and recent environmental changes have made a positive contribution to expanding greenery in Nepal.In general, NDVI has increased in more than 80% of the study area with a decreasing trend in some of the pixels in both southern and northern parts.In the south, the population density and cultivated land are high, but in the north, the area is mostly covered by ice and glaciers that result in low NDVI.We also examined the NDVI variations in different altitudinal zones.The hills have high NDVI values due to the presence of a large number of community forestry The forest biomass carbon density showed spatial heterogeneity (Figure 6b).It was high in the south-western and central areas but lower in eastern and northern parts of Nepal.The highest carbon density reached 168.82 t C/ha in western parts and the lowest was 9.27 t C/ha in upper parts of hills mainly occupied by sub-alpine forests.In Figure 6b, the green color represents the dense broad leaved forest.The light green colors also represent the broad leaved open and closed forest with few needle leaved forest.Similarly, the red and yellow colors mostly represent needle leaved forests in Nepal (Figure 1).

Discussion
The NDVI experienced increasing trends in large areas of Nepal from 1982-2015 (Table 1; Figures 2  and 4).The annual NDVI increased by 0.0008 year −1 with positive trends in the estimation of both the linear and Sen's slope for all seasons, with obvious spatial heterogeneity during the last three decades.The majority of the area has positive NDVI trends with negative trends in a few locations.The results are consistent with previous studies carried out in several parts of the world, mainly in the Northern Hemisphere [17,19,21,78,84].Vegetation change is regarded as an indicator of environmental change which has severe ecological consequences in Nepal.The ecology of vegetation changes has significance impacts on functioning ecosystems, energy balance, species richness and land atmospheric interactions.Generally, high NDVI values (indicated by an integrated NDVI or the annual maximum NDVI value) reflect a good status of plant productivity and biomass.A positive NDVI trend indicates that ecology and ecosystem conditions are better.The ecological restoration and recent environmental changes have made a positive contribution to expanding greenery in Nepal.In general, NDVI has increased in more than 80% of the study area with a decreasing trend in some of the pixels in both southern and northern parts.In the south, the population density and cultivated land are high, but in the north, the area is mostly covered by ice and glaciers that result in low NDVI.We also examined the NDVI variations in different altitudinal zones.The hills have high NDVI values due to the presence of a large number of community forestry and agricultural practices.The tree line ecotone ranges between 3000-4300 m altitude [85] and also has a good NDVI distribution.It can be said that the tree line area is experiencing good ecological restoration with more seedlings and saplings [86].The NDVI distribution and trends in different ecological zones indicate a good ecological status with a positive level of photosynthetic activities, biomass, productivity and ecosystem balance [19,23,87,88], as well as a lengthening of the growing season [89][90][91].
The analysis of the climate data showed that annual average temperature has significantly increased but the annual average precipitation has decreased.In our examination, NDVI and temperature have a positive correlation at annual and seasonal time scales except in the pre-monsoon season in which the correlation is negative.The positive correlation between NDVI and temperature across a large area with a high fraction of significant correlation indicates that temperature is one of the major factors driving vegetation growth.Conversely, the annual and seasonal correlation between NDVI and precipitation is negative, except for a positive correlation in winter and post-monsoon seasons.The positive and significant correlation between NDVI and climatic controls as shown in Figure 5 showed that climate is responsible for increasing vegetation: 34.97% of the increased vegetation is closely related to temperature (R T is significant in 34.97% of the area) and 1.08% of the increased vegetation is related to precipitation (R p is significant in 1.08% area).The correlation between NDVI and precipitation was mostly negative at seasonal and annual time scales, but the correlation in each pixel was positive in 30.76% of the area, with 1.08% being significant and 16.35% being significant negatively.The correlation between NDVI and precipitation has a lower significantly positive and higher significantly negative correlation than the correlation between NDVI and temperature in Nepal.In Nepal, a change in NDVI was positively correlated with temperature in the majority of the area compared with precipitation.Mostly, in the mid-hills, the correlation between NDVI and temperature was significant and there was a significant negative correlation between NDVI and precipitation in south-western parts of Nepal.The western parts of Nepal are located far from the ocean relative to the eastern parts, so that the rainfall was less but the vegetation distribution and presence of forest carbon were better than in eastern Nepal.Similar types of relationships, i.e., positive with respect to temperature and negative with respect to precipitation between NDVI and climate were reported in global and regional-scale research [16,22,84,92].Besides the temperature and precipitation, the vegetation can also be influenced by other factors in Nepal which are related to CO 2 [44], topography and human disturbance [85], moisture limitation [93], nitrogen deposition, irrigation, change in croplands, land use practices and ecological restoration.A previous study discovered positive NDVI trends in Nepal which were mainly attributed to CO 2 fertilization [44].In our study, we also found positive NDVI trends and a strong correlation with temperature compared to the precipitation.In addition, we derived positive BCD (t C/ha) trends based on NDVI which is added value and different from the study carried out by [44] in Nepal.Nepal is a country of diverse topography, with a high climatic gradient and complex network of ecosystems where multiple factors might influence vegetation changes.We have recognized the limitation of our study that focuses on climatic variables as potential driving factors because the temperature and precipitation change rapidly in Nepal.The study of vegetation changes associated with several environmental factors indicates that vegetation changes may be related to multiple factors which require careful investigation with high quality input in future.
Various researchers developed different empirical models to estimate global forest biomass using satellite-derived NDVI and geographical location [30,32,94,95].Among them, we followed the model used by [32] because this model used NDVI max to predict the biomass C density, which is also suitable for Nepal, rather than using growing season and annual average NDVI.In addition, this model has a high prediction accuracy indicated by a higher R 2 of 0.64 (p < 0.001) compared to the model using growing season and annual average NDVI (R 2 = 0.53, p < 0.001).We used this model as Nepal belongs to the same hemisphere and similar regions near to Himalaya assuming similar patterns of forests, topography, and climates.Moreover, the result obtained was similar to the government's carbon stock inventory of Nepal.Our estimation found an average of 115.392 t C/ha but the government's field-based study obtained an average value of 108.62 t C/ha from 2010-2014 [53].The present study indicates an increase of 6.77 t C/ha in average carbon stock difference of 6.23% compared to that of the government's inventory data.The forest inventory data considered the trees above 10 cm diameter at breast height (DBH) only, whereas RS-based carbon estimation considers vegetation of all sizes, but the RS-based method may result in slightly higher values of carbon stock than those of a field-based inventory.The present carbon study is limited to forest land only, which occupied 40.36% of the total area of Nepal [53], as shown in Figure 1.The RS-based carbon estimation helps to deal with the issues raised with respect to the reduction of greenhouse gas emissions committed to during the Kyoto Protocol United Nations Framework Convention on Climate Change (UNFCC) and ecosystem functioning in Nepal.
The validity and accuracy of the results were checked based on a comparison with several global, regional, plateau and national-scale research and government reports of Nepal.As mentioned above, the result is more consistent with previous research carried out in the Northern Hemisphere, China [1,84,96] and Tibetan plateau [21], where the NDVI is reportedly positive and mostly related to temperature and precipitation.The results were logical with respect to different land use types [62], vegetation types [53] and climatic variability [97] in Nepal.The tropical zone (<1000 m altitude) in eastern Nepal is covered with evergreen forests but in the tropical zone of western Nepal, there is dense forest.In the temperate and sub-alpine zone (2000-4000 m) of Nepal, there are moist temperate deciduous forests and evergreen coniferous forests from west to east.The alpine zones (above 4000 m altitude) are dominated by moist and dry alpine scrub above the timber line [53].The alpine meadows which are mainly grazing lands are found in alpine areas, especially in the trans-Himalayan zone.These different vegetation types in different regions lead to NDVI differences in Nepal (Figures 2 and 4).The afforestation program after the initiation of the community forestry project in 1980 has changed the regional vegetation composition.Before the 1980s, it was more forest and grasslands in the western parts of Nepal than in the eastern and central parts of Nepal; eastern Nepal has the most crop lands followed by central Nepal.After three decades, the NDVI distribution and changing patterns are spatially and temporally varied (Figures 2 and 4).Similarly, the obtained results of the forest carbon biomass were also verified with national forest inventory reports of Nepal [53,98].In this study, we used a proxy NDVI for the analysis.It has a low resolution but the longest temporal coverage and the best data set available for long-term analysis of vegetation dynamics [59].These data sets have been widely used in the studies of vegetation changes because of the high quality and elimination of possible errors caused by integral and extraneous factors [2,57,58,99].This data set has also proved to be one of the best products to detect temporal changes in vegetation [60].This research is based on country-wide data at large temporal and spatial scales.The leaf area index, solar radiation, soil moisture and other physical and eco-climatological factors were not used in this analysis.The use of ground-based ecological data and ecosystem modeling could also be a good analysis in further research.

Conclusions
The study provides spatial and temporal NDVI trends and the correlation between NDVI and climatic factors on annual and seasonal time scales during the last three decades in Nepal.It further provides forest BCD (t C/ha) and its temporal trends in Nepal.We found positive NDVI and carbon trends during the study period in which annual, pre-monsoon and monsoon season NDVI experienced a significant positive trend while winter and post-monsoon NDVI experienced positive but insignificant trends.Spatially, the majority of the areas experienced significant positive NDVI trends in Nepal.The NDVI distribution and positive NDVI changing ratio were higher in the hills and middle parts of Nepal.The NDVI-based carbon analysis showed increasing trends in forests with higher biomass carbon density in western parts than in eastern parts of the mid-hills.We found that increased NDVI was mostly correlated with increasing temperature rather than precipitation.Thus, temperature is recognized as a principal climatic factor affecting NDVI but not limited to Nepal.The overall findings of this study aid in understanding the relationship between climate and long-term NDVI changes, and the implication of NDVI for carbon dynamics helps to provide a scientific knowledge for ecology and environmental management in Nepal.

Figure 1 .
Figure 1.Study area: (a) Location map; (b) Ecological zone and administrative boundary of the study area; (c) land use types in Nepal.

Figure 1 .
Figure 1.Study area: (a) Location map; (b) Ecological zone and administrative boundary of the study area; (c) land use types in Nepal.
, and the annual Kendall tau of NDVI was 0.45 from 1982-2015.The positive Kendall's tau statistics also indicate a positive NDVI trend.The figures above were developed based on obtained Mann-Kendall test statistics and Sen's slope values which showed annual and seasonal NDVI trends and inter-annual variation of mean NDVI distribution in Nepal from 1982 to 2015.The analysis of temporal NDVI changes in Terai (<1000 m), Hill (1000-3000 m) and Mountain (>3000 m) areas showed that NDVI has highest in the hills, with an annual positive trend of 0.0012, followed by Terai with a mean NDVI of 0.56 and annual average trend of 0.001 (Table 4 to 0.5 in the low lands of southern regions.The maximum mean NDVI ranges from 0.7 to 0.8, mainly in mid-hill areas occupied by broadleaved dense forest.The average NDVI values in western Nepal (from the western boundary to 83 • E) and central Nepal (from 83 • to 86 • 3 E) were 0.481 and 0.485, respectively, but 0.551 in eastern Nepal (from 86 • 3 E to eastern boundary).

Figure 4 .
Figure 4. NDVI trends and significance of the NDVI trends at annual and seasonal time scales: (a) Annual NDVI; (b) winter season NDVI; (c) pre-monsoon season NDVI; (d) monsoon season NDVI and (e) post-monsoon season NDVI trend (year −1 ) in Nepal, 1982-2015 (small maps at the bottom of annual and season slope maps show p-values of NDVI trends).

Figure 4 .
Figure 4. NDVI trends and significance of the NDVI trends at annual and seasonal time scales: (a) Annual NDVI; (b) winter season NDVI; (c) pre-monsoon season NDVI; (d) monsoon season NDVI and (e) post-monsoon season NDVI trend (year −1 ) in Nepal, 1982-2015 (small maps at the bottom of annual and season slope maps show p-values of NDVI trends).

Figure 5 .
Figure 5. Spatial and temporal correlation between: (a) NDVI and temperature (R T ); (b) NDVI and precipitation (R P ); (c) significant correlation between NDVI and temperature; (d) significant correlation between NDVI and precipitation in Nepal, 1982-2015.

Figure 6 .
Figure 6.The temporal and spatial variation of BCD (t C/ha) in Nepal: (a) temporal variation of BCD (t C/ha); (b) spatial variation of forest BCD (t C/ha) in Nepal, 1982-2015.

Figure 6 .
Figure 6.The temporal and spatial variation of BCD (t C/ha) in Nepal: (a) temporal variation of BCD (t C/ha); (b) spatial variation of forest BCD (t C/ha) in Nepal, 1982-2015.

Table 2 .
Annual and seasonal mean NDVI and trends (year −1 ) in three different ecological zones in Nepal, 1982-2015.

Table 2 .
Annual and seasonal mean NDVI and trends (year −1 ) in three different ecological zones in Nepal, 1982-2015.

Table 3 .
Annual and seasonal correlation coefficients between NDVI and temperature (R T ) and precipitation (R P ) in Nepal from 1982-2015; the values in parentheses refer to the p values of the correlation coefficients.