An Assessment of Environmental Impacts on the Ecosystem Services: Study on the Bagmati Basin of Nepal

: The upsurges in population, internal migration, and various development works have caused signiﬁcant land use and land cover (LULC) changes in the Bagmati Basin of Nepal. The e ﬀ ects of climate change such as increased precipitation and temperature are a ﬀ ecting the provision of ecosystem services (ES). In this regard, this study particularly treated water yield (WY), soil loss, nitrogen export, and carbon ﬂuctuation in the basin. Integrated Valuation of Ecosystem Services and Tradeo ﬀ s (InVEST) tools were used to carry out a comparative analysis of ES based on LULC data for 2000 and 2010 and corresponding climate data. To analyze the future period (2010–2099), we have used climate data from the multi-model ensemble (MME) of statistically downscaled and bias-corrected 12 best global climate models (GCMs). A raw GCM analysis (based on historical observational data) from 29 GCMs was done ﬁrst. The results shows with a subsequent degradation of ES providers like forests and an increment in agricultural and urban areas, ES are on a verge of degradation. Furthermore, a projection of future climate patterns depicts increased precipitation and temperature. Thus, urgent measures are required for the sustainable provision of ES. Outcomes of the study are expected to help in the incorporation of ES in development policies promoting low-impact development along with maintaining ecological and economic goals. The study closes by presenting a recommendation for model application and future study needs.


Introduction
Ecosystem services (ES) are defined as services provided by nature that directly and indirectly benefit the well-being of humans. The concept of ES has been developed in the scientific literature since the end of the 1970s [1], and the Millennium Ecosystem Assessment (MEA) report 2005 [2] is considered the milestone in the field of mapping ES. The MEA report evaluated the impacts on human beings as the consequences of changes in the ecosystem. The core finding of the MEA report is that human activities are depleting the Earth's natural resources, which is inducing stress on the self-sustaining capability of the natural environment. This phenomenon urges immediate measures for checking and prevention in order to enhance the sustainability of ecosystems for future generations.
To meet the requirements of a globally increasing population for food, shelter, timber, etc., land use changes are occurring but at the expense of a degraded environment and ecosystems. Land use and land cover (LULC) change is considered the major factor for the degradation of ES and loss in biodiversity [3,4]. Climate change is another factor impacting natural systems, and ultimately affects the flow of ES [5]. With the increasingly evident effects of climate change, which is especially impacting

Study Area
The Bagmati Basin lies in the middle mountain region of Nepal ( Figure 1) at latitude 26°42′ to 27°50′ N and longitude 85°2′ to 85°58′ E. The Bagmati River is the principal river of the basin and it originates from the north of the capital city, Kathmandu, at Shivpuri (Bagdwar) at an altitude of 2690 m. The river flows through Kathmandu and provides most of the city's drinking water in the upper part of the basin, supports hydropower generation in the middle basin and in the lower part of the basin, and supports large-scale irrigation for agriculture [29]. The climate of the basin varies with elevation; higher mountains have a cold temperate climate, mid-elevation levels have a warm temperate climate, and the southern lowlands have a subtropical climate [30]. The basin has a mean annual temperature of 20-30 °C and the mean annual precipitation is about 1800 mm [30]. The Bagmati River is a spring and monsoon rain-fed river and consists of many tributaries such as the Hanumante, Manohara, Dhobi Khola, etc. The total area of the basin is 375,000 ha and, owing to the good data availability for this study, a 276,897 ha area was selected. This represents the area above the Pandredhovan gauge station ( Figure 2). We chose the Bagmati Basin as our study area as it is facing huge anthropogenic pressures associated mostly with centralized facilities, urban expansion, and increased agricultural activity.   Economic Cooperation Climate Centre (APCC) was used for statistical downscaling of bias-corrected global climate model (GCM) data.

Study Area
The Bagmati Basin lies in the middle mountain region of Nepal ( Figure 1) at latitude 26°42′ to 27°50′ N and longitude 85°2′ to 85°58′ E. The Bagmati River is the principal river of the basin and it originates from the north of the capital city, Kathmandu, at Shivpuri (Bagdwar) at an altitude of 2690 m. The river flows through Kathmandu and provides most of the city's drinking water in the upper part of the basin, supports hydropower generation in the middle basin and in the lower part of the basin, and supports large-scale irrigation for agriculture [29]. The climate of the basin varies with elevation; higher mountains have a cold temperate climate, mid-elevation levels have a warm temperate climate, and the southern lowlands have a subtropical climate [30]. The basin has a mean annual temperature of 20-30 °C and the mean annual precipitation is about 1800 mm [30]. The Bagmati River is a spring and monsoon rain-fed river and consists of many tributaries such as the Hanumante, Manohara, Dhobi Khola, etc. The total area of the basin is 375,000 ha and, owing to the good data availability for this study, a 276,897 ha area was selected. This represents the area above the Pandredhovan gauge station ( Figure 2). We chose the Bagmati Basin as our study area as it is facing huge anthropogenic pressures associated mostly with centralized facilities, urban expansion, and increased agricultural activity.

Model Description and Data
For the assessment of ES provision in the basin, four ES, namely, WY, soil loss, nitrogen export, and carbon storage, were evaluated on the sub-basin scale in the Bagmati Basin of Nepal. We have done a comparative study of ES for the years 2000 and 2010 based on the LULCs of 2000 and 2010 and corresponding climate data and for a future period, ES were mapped based on the LULC map of 2010 and projected climate data. The LULC map was obtained from the ICIMOD Nepal Geospatial Portal (http://rds.icimod.org/) which is prepared using public domain Landsat Thematic mapper(TM) data. Uddin et al. [31] have discussed the detailed description of the data including LULC classification, accuracy assessment, and limitations, hence we have used the LULC data from the ICIMOD portal referring to the findings from Uddin et al. [31]. Likewise, the rainfall and temperature data for the period of 1987-2016 from 25 stations from all over Nepal (Table S1, Supplementary Materials) have been used to downscale GCM data for Nepal. The missing data were extracted from MERRA2 grid data [32] after doing quality control (checking daily correlation coefficient) using R package. GCM data from 10 stations (stations in bold in Table S1, Supplementary Materials) have been used for a future period (2010-2099) ES evaluation of the Bagmati Basin. As well, for the analysis of the Bagmati Basin for the period 2000 and 2010, a precipitation raster is prepared using data from 23 rain stations ( Figure 3) and, for preparation of an evapotranspiration raster, precipitation and temperature data from 9 stations (stations in bold in Table S1

Model Description and Data
For the assessment of ES provision in the basin, four ES, namely, WY, soil loss, nitrogen export, and carbon storage, were evaluated on the sub-basin scale in the Bagmati Basin of Nepal. We have done a comparative study of ES for the years 2000 and 2010 based on the LULCs of 2000 and 2010 and corresponding climate data and for a future period, ES were mapped based on the LULC map of 2010 and projected climate data. The LULC map was obtained from the ICIMOD Nepal Geospatial Portal (http://rds.icimod.org/) which is prepared using public domain Landsat Thematic mapper(TM) data. Uddin et al. [31] have discussed the detailed description of the data including LULC classification, accuracy assessment, and limitations, hence we have used the LULC data from the ICIMOD portal referring to the findings from Uddin et al. [31]. Likewise, the rainfall and temperature data for the period of 1987-2016 from 25 stations from all over Nepal (Table S1, Supplementary Materials) have been used to downscale GCM data for Nepal. The missing data were extracted from MERRA2 grid data [32] after doing quality control (checking daily correlation coefficient) using R package. GCM data from 10 stations (stations in bold in Table S1, Supplementary Materials) have been used for a future period (2010-2099) ES evaluation of the Bagmati Basin. As well, for the analysis of the Bagmati Basin for the period 2000 and 2010, a precipitation raster is prepared using data from 23 rain stations ( Figure 3) and, for preparation of an evapotranspiration raster, precipitation and temperature data from 9 stations (stations in bold in Table S1  For the downscaling of GCMs for future climate data, we have used APCC's AIMS software. We have used the multi-model ensemble (MME) of 12 best GCMs of Coupled Model Intercomparison Project 5 (CMIP5) under two Representative Concentration Pathways (RCP) scenarios namely RCP 4.5 and RCP 8.5 to study the projection of ES in the future. The descriptions of RCP 4.5 and RCP 8.5 are given in Figure 4. InVEST tools were used to map WY, nitrogen export, and carbon storage, and a RUSLE method was used to map soil loss. As carbon storage computation by the InVEST carbon model is independent of the climate data, carbon storage was evaluated with only LULC change and the corresponding climate scenario whereas WY, soil loss, and nitrogen export were evaluated for both scenarios of LULC change and future climate change scenarios. We choose the InVEST model in our study due to its capability to give a quick overview of ES, especially in regions with low data availability [33] in developing nations like Nepal. InVEST is a suite of a free and open source software models. Its modular design provides an effective tool for evaluating the possible outcomes of For the downscaling of GCMs for future climate data, we have used APCC's AIMS software. We have used the multi-model ensemble (MME) of 12 best GCMs of Coupled Model Intercomparison Project 5 (CMIP5) under two Representative Concentration Pathways (RCP) scenarios namely RCP 4.5 and RCP 8.5 to study the projection of ES in the future. The descriptions of RCP 4.5 and RCP 8.5 are given in Figure 4. InVEST tools were used to map WY, nitrogen export, and carbon storage, and a RUSLE method was used to map soil loss. As carbon storage computation by the InVEST carbon model is independent of the climate data, carbon storage was evaluated with only LULC change and the corresponding climate scenario whereas WY, soil loss, and nitrogen export were evaluated for both scenarios of LULC change and future climate change scenarios. We choose the InVEST model in our study due to its capability to give a quick overview of ES, especially in regions with low data availability [33] in developing nations like Nepal. InVEST is a suite of a free and open source software models. Its modular design provides an effective tool for evaluating the possible outcomes of alternative management and climate scenarios and it helps in the selection of the best decision scenarios [27].
Sustainability 2020, 12, x FOR PEER REVIEW 5 of 22 alternative management and climate scenarios and it helps in the selection of the best decision scenarios [27].

Water Yield Estimation
The InVEST water yield model was used for the WY estimation and this model is based on the water balance principle and Budyko curve. It determines the amount of water running off each pixel as the precipitation minus the fraction of the water that undergoes evapotranspiration. The input data included precipitation, evapotranspiration, LULC, soil properties, basin boundary, and biophysical attributes ( Table 1). The spatial resolution of LULC, soil depth, and plant-available water content used in the study was 30 m. The annual WY Y(x) on each pixel of the landscape (x) is determined as follows: AET(x) = Annual actual evapotranspiration for pixel x P(x) = Annual precipitation on pixel x.

Water Yield Estimation
The InVEST water yield model was used for the WY estimation and this model is based on the water balance principle and Budyko curve. It determines the amount of water running off each pixel as the precipitation minus the fraction of the water that undergoes evapotranspiration. The input data included precipitation, evapotranspiration, LULC, soil properties, basin boundary, and biophysical attributes ( Table 1). The spatial resolution of LULC, soil depth, and plant-available water content used in the study was 30 m. The annual WY Y(x) on each pixel of the landscape (x) is determined as follows: AET(x) = Annual actual evapotranspiration for pixel x P(x) = Annual precipitation on pixel x.
AET(x) P(x) is the evapotranspiration portion of the water balance for vegetated LULC types. This is computed based on an expression of the Budyko curve proposed by Fu [34] and Zhang et al. [35]: where PET(x) = potential evapotranspiration calculated using Equation (3); ω(x) = non-physical parameter (characterizes the natural climatic-soil properties). Calculation of PET is as follows: where ET o (x) = reference evapotranspiration from pixel x. Kc(l x ) = plant evapotranspiration coefficient, which is associated with the LULC l x on pixel x. ω(x) is an empirical parameter. It is determined as the expression given by Donohue et al. [36]: where ω(x) is based on plant-available water content (PAWC), precipitation, and the Z constant. The Z constant defines the local precipitation pattern and additional hydrogeological characteristic of the basin [27].
The temperature-based Hargreaves equation was used for the computation of reference evapotranspiration, given that this method generates superior results to the Pennman-Montieth method in the case when long-term data is limited [37]. An average soil depth and PAWC raster map was prepared using the Soil and Terrain (SOTER) database of Nepal and ArcGIS. The biophysical information on the LULC code and its descriptive names, the maximum root depth for vegetated land use classes in mm, and the plant evapotranspiration coefficient for each LULC class was required ( Table 2). In this study, the root depth of main vegetation types was obtained following Chen et al. [38] and the evapotranspiration coefficient of each LULC type used on the model was based on Allen et al. [39] and the InVEST user manual [27].

Soil Loss Estimation
The RUSLE model coupled with ArcGIS was used for the computation of soil loss in the basin. Due to data simplicity and the provision of scenario analysis and taking measures against erosion, RUSLE is widely used at large scales for soil loss assessment [40]. It computes soil erosion as the product of six factors representing rainfall erosivity, soil erodibility, slope length, slope steepness, cover and management practices, and supporting conservation practices: where A = average annual soil loss amount in (t/ha/yr) R = rainfall-runoff erosivity factor (MJ mm/h/ha/yr) K = soil erodibility factor L = slope length factor S = slope steepness factor C = land cover management factor P = support practice factor Among several equations available for the rainfall-runoff erosivity factor (R), we have used the equation by Singh et al. [41] as it is the one recommended for the Himalayan region: where R = rainfall erosivity factor (MJ mm ha −1 h −1 yr −1 ) P = mean annual precipitation (mm).
The K factor is determined as per the Food and Agriculture Organization (FAO) database adapted to Nepal by World Soil Information (ISRIC). The soil unit map was extracted from the SOTER database of Nepal and the K factor was computed based on different published literature on mountainous areas [42,43] and other countries [44]. The K factor of the study area varied from 0.19 to 0.49.
The slope formula based on slope length was used for computation of slope length factor (L) based on references [45,46] which is given as: where λ indicates a field slope length and 22.13 is the slope length of a unit runoff plot (m). m = slope length exponent. The slope steepness factor (S) represents the effect of slope steepness on the intensity of soil erosion. The factor is calculated using Equation (8) as described by Wischmeier and Smith [45]: where s = slope in percent, which is determined from a digital elevation model (DEM). Likewise, the cover management factor (C) and the support practice factor (P) were obtained from various published reports [25,47] in a similar area. Separate raster layers for each factor were prepared in GIS and the average annual soil erosion rate was determined by multiplying the respective factors in the ArcGIS environment.

Nitrogen Export Estimation
The nutrient delivery ratio (NDR) model was used to map nitrogen export in the basin. The model uses a simple mass balance approach and describes the movement of nutrient mass through space and aims to quantify nutrient (nitrogen and phosphorous) export. The model maps the flow of nutrients from various sources to the stream network. Sources of nutrients are determined based on a LULC map and associated loading rates. The nutrient loads are divided into sediment-bound and dissolved parts, which will be carried by surface and subsurface flow respectively and stops when they reach a stream. Next, the delivery factors were computed for each pixel, based on the properties of pixels belonging to the same flow path (their slope and retention efficiency of each land use) [27]. All pixel-level contributions were summed at the basin/sub-basin outlet to compute the nutrient export: where X expton = total export amount of nutrients in the basin (kg. yr −1 ) X expi = export amount of nutrients from each grid NDR Subs, i = subsurface nutrient transfer rate.
The raster layers required as an input for running the model are raster maps of the DEM, LULC, and precipitation. The model requires the biophysical table having a coefficient for nitrogen loading for each LULC category ( Table 2). These values were obtained following Sharp et al. [27] and Line et al. [48].

Carbon Storage Estimation
The InVEST carbon model used in this study maps carbon storage densities to the LULC raster by aggregating the amount of carbon stored on four major carbon pools. The carbon pools considered were aboveground biomass, belowground biomass, soil, and dead organic matter to produce the total amount of carbon storage.
The carbon storage C m,i,j for a given grid cell (i, j) with land use type m can be calculated as: where A = actual area of each grid cell (ha); Ca m,I,j , Cb m,I,j , Cs m,I,j , Cd m,I,j are the aboveground, belowground, soil organic, and dead organic matter carbon density in MgC/ha for grid cell (i, j) with land use type m ( Table 2). The carbon pool data were obtained from published literature [49,50] and InVEST user guidelines.

GCM Downscaling and Future Climate Scenarios
For this study, 29 GCMs of CMIP5 were statistically downscaled and bias-corrected using the simple quantile mapping (SQM) [51] method. We adopted RCP 4.5 and RCP 8.5 concentration pathways (Table 3, adapted from reference [52]) for this study. The future projections of daily precipitation and temperature were performed at 25 meteorological stations of Nepal. Statistical downscaling makes it possible to perform a quantitative comparison with observational data through bias correction using the observations and it is easy to convert them into high-resolution data. The multi-model ensemble was prepared for the study basin using 12 best GCMs after doing raw GCM analysis of the downscaled GCMs. For the downscaling, we used the AIMS developed by the APCC. The AIMS module is a free and open source module available online from www.aims.apcc21.org [53].

Evaluation of Land Use and Land Cover Change (LULC)
A total of 37,487.32 ha (13.54% of the study area) of land cover has faced conversion from one land use type to another in the period between 2000 and 2010. Among all the conversions from one type to another, conversion to agriculture and built-up areas is highest on the basin (Tables 4 and 5). The highest overall increment is in the agricultural area (increased by 9232 ha), followed by the built-up area (increased by 4087.28 ha). Similarly, the highest reduction is on grassland which decreased by 6012.12 ha. The built-up area in 2010 had increased by 4087.28 ha compared to 2000 and in this conversion, 82.71% (3380.96 ha) came from the agricultural area only. This decrement in agricultural area is mostly concentrated around major cities of the Kathmandu valley, changing the major traditional agricultural pattern in the valley. Likewise, 9362.92 ha of forest land and 7555.8 ha of grassland were converted to agricultural areas from 2000 to 2010. However, as the conversion from agricultural to other land use types like built-up and forest areas is notable, the net gain on agricultural areas was 9232.00 ha in 2010. The grassland land use faced the highest decrement contributing to other land use types mostly for the agriculture land use (7555.8 ha) and barren land has also faced similar conversion with the highest (3194.24 ha) conversion to agricultural land use types. The upper part of the basin has faced significant expansion on urban and agricultural areas and this transition has resulted due to several factors like internal migration of people from other parts of the country, increased agricultural practices, ongoing infrastructural activities, and economic movements in the basin.

Climate Data Analysis
Twelve GCMs (Table 6) were selected for the study basin after raw GCM analysis ( Figure 5) compared to the historical observation data (data period 1987-2016). With reference to the findings from several studies, multi-model ensemble (MME) accounted for the uncertainty of the single model [51,54], MME was prepared with an ensemble average of 12 best GCMs for the study area. The data were analyzed based on six scenarios, three each for RCP 4.5 and RCP 8.5 (Table 7).     The patterns of two climate parameters (temperature and precipitation) were observed for the period 2010-2099 compared with the baseline historical period of 1987-2016 ( Table 7). The average annual precipitation for the historical period was 1609.57 mm. It was observed that the average annual precipitation is expected to increase by 94.

Water Yield
The average annual precipitation was 1879.53 mm and 1766.59 mm for the periods of 2000 and 2010, respectively. Comparing the two periods, the WY is observed to decrease by 106.592 MCM (million cubic meters) from 2000 to 2010. Figure 8 shows a raster map for the difference of WY between 2000 and 2010. In 2000, the total average annual WY was observed to be 3278.609 MCM and in 2010, WY was 3172.017 MCM. Performing a sub-basin-scale evaluation, WY was observed to have decreased in 2010 in all sub-basins except in sub-basin 2 ( Table 8). As the computation of WY is a factor of evapotranspiration along with precipitation and LULC change, increment in the water yield in sub-basin 2 is attributable to an increase in the built-up area, which comes up with several other consequences like increased nitrogen load and reduced carbon storage.
scenarios. The WY per ha is highest on sub-basin 5 and lowest on sub-basin 1 in 2010 but, in future scenarios, sub-basin 6 will have the highest WY and sub-basin 2 will have the lowest water yield per ha of land (Table 9). This provides a general estimate and trend of WY but as WY is a function of evapotranspiration, highly affected by LULC, projections based on future LULC scenarios can give a better estimate.   Table 9. Water yield per ha computation for period S1-S3 under RCP 4.5 and RCP8.5 (m 3 /ha).  The assessment of WY on climate scenarios based on the baseline land use of 2010 and the corresponding climate data obtained from MME of GCMs shows increased WY with an increased projection of precipitation in future scenarios (Table 9). Compared to 2010, the total WY of the basin is observed to decrease in S1Rcp4.5 and S1Rcp8.5. However, even with the baseline LULC of 2010, the total WY was observed to increase from period S1 to S2 to S3 under both RCP 4.5 and RCP 8.5 scenarios. The WY per ha is highest on sub-basin 5 and lowest on sub-basin 1 in 2010 but, in future scenarios, sub-basin 6 will have the highest WY and sub-basin 2 will have the lowest water yield per ha of land (Table 9). This provides a general estimate and trend of WY but as WY is a function of evapotranspiration, highly affected by LULC, projections based on future LULC scenarios can give a better estimate.

Soil Loss
The annual average soil loss is reduced in 2010 (20.46 Mt/yr) compared to 2000 (21.38 Mt/yr) and is attributable to decreased annual average precipitation. Figure 9 shows a raster map for the difference of soil loss between 2000 and 2010. In 2000, agriculture accounted for the highest average rate of soil loss (198.40 t/yr/ha) whereas, in 2010, barren land was responsible for the highest soil loss (225.62 t/yr/ha) ( Table 10). In both periods, the highest average rate of soil loss was as a result of agricultural land use causing 56.07% and 56.06% of total soil loss in 2000 and 2010, respectively.

Sub-Basin
2010 S1_Rcp4.5 S2_Rcp4.5 S3_Rcp4.5 S1_Rcp8.5 S2_Rcp8.5 S3_Rcp8.5 38 Mt/yr) and is attributable to decreased annual average precipitation. Figure 9 shows a raster map for the difference of soil loss between 2000 and 2010. In 2000, agriculture accounted for the highest average rate of soil loss (198.40 t/yr/ha) whereas, in 2010, barren land was responsible for the highest soil loss (225.62 t/yr/ha) ( Table 10). In both periods, the highest average rate of soil loss was as a result of agricultural land use causing 56.07% and 56.06% of total soil loss in 2000 and 2010, respectively.
The assessment of soil loss on future periods under RCP 4.5 and RCP 8.5 was computed with the baseline LULC map of 2010 and projected MME of climate data. The increment in precipitation will result in increased rainfall erosivity and, hence, soil loss is projected to increase in all consequent future scenarios. Soil loss increased by 1.53, 1.74, and 3.1 Mt/yr in S1, S2, and S3 periods, respectively, under the RCP 4.5 scenario. Likewise, under RCP 8.5, soil loss was projected to increase by 3.63, 4.23, and 6.59 Mt/yr in S1, S2, and S3 periods, respectively (Table 11).    The assessment of soil loss on future periods under RCP 4.5 and RCP 8.5 was computed with the baseline LULC map of 2010 and projected MME of climate data. The increment in precipitation will result in increased rainfall erosivity and, hence, soil loss is projected to increase in all consequent future scenarios. Soil loss increased by 1.53, 1.74, and 3.1 Mt/yr in S1, S2, and S3 periods, respectively, under the RCP 4.5 scenario. Likewise, under RCP 8.5, soil loss was projected to increase by 3.63, 4.23, and 6.59 Mt/yr in S1, S2, and S3 periods, respectively (Table 11).

Nitrogen Export
Nitrogen export from the various LULC types directly affects the watershed health, humans, and aquatic life processes. The comparative study (Table 12) of nitrogen export between LULC of 2000 and 2010 ( Figure 10) showed that the total nitrogen export of the basin increased by 34,490.90 kg in 2010 compared to 2000. The sub-basin-scale evaluation shows that nitrogen export was the highest from sub-basin 2 in both periods. The assessment of nitrogen export on future scenarios with projected MME of climate data and the baseline LULC of 2010 did not show a precise pattern like other ES. The nitrogen export was determined based on LULC and associated loading rates. As the study was conducted on baseline land use of 2010, the comparative value of nitrogen export on climate scenarios was observed to be less than that of the 2010 case. The estimation with future land use scenarios can enhance the accuracy of estimation. However, in all cases (Table 13), nitrogen export was highest in sub-basin 2 and lowest in sub-basin 7 like in the baseline period of 2010. highest from sub-basin 2 in both periods. The assessment of nitrogen export on future scenarios with projected MME of climate data and the baseline LULC of 2010 did not show a precise pattern like other ES. The nitrogen export was determined based on LULC and associated loading rates. As the study was conducted on baseline land use of 2010, the comparative value of nitrogen export on climate scenarios was observed to be less than that of the 2010 case. The estimation with future land use scenarios can enhance the accuracy of estimation. However, in all cases (Table 13), nitrogen export was highest in sub-basin 2 and lowest in sub-basin 7 like in the baseline period of 2010.

Carbon Storage
The total modeled carbon was reduced by 969,923.33 Mg from 2000 to 2010. The reduction in forest cover, shrubland, and grassland caused the overall reduction of carbon storage on the basin. Likewise, a sub-basin-scale evaluation of carbon storage shows that, in both periods, carbon storage per ha was the highest in sub-basin 7 and lowest in sub-basin 2 (Table 14). Figure 11 shows a raster file for the difference in carbon storage between 2000 and 2010. There were increased anthropogenic activities in sub-basin 2 which is mostly comprised of residential/built-up areas and agricultural areas of the basin. This represents the area with the lowest storage despite having the highest area. Likewise, comparing total carbon storage in 2000 and 2010, sub-basin 2 marks the highest loss with a 315,381.00 Mg reduction and sub-basin 6 has increased storage by 168,886.91 Mg.

Carbon Storage
The total modeled carbon was reduced by 969,923.33 Mg from 2000 to 2010. The reduction in forest cover, shrubland, and grassland caused the overall reduction of carbon storage on the basin. Likewise, a sub-basin-scale evaluation of carbon storage shows that, in both periods, carbon storage per ha was the highest in sub-basin 7 and lowest in sub-basin 2 (Table 14). Figure 11 shows a raster file for the difference in carbon storage between 2000 and 2010. There were increased anthropogenic activities in sub-basin 2 which is mostly comprised of residential/built-up areas and agricultural areas of the basin. This represents the area with the lowest storage despite having the highest area. Likewise, comparing total carbon storage in 2000 and 2010, sub-basin 2 marks the highest loss with a 315,381.00 Mg reduction and sub-basin 6 has increased storage by 168,886.91 Mg.

Discussion
The global concern for the conservation and promotion of ES is rising along with an awareness that projected extreme events as a result of climate change will hamper the provision of services. Normally, ES are renewable if they can be managed sustainably but can be depleted or degraded if mismanaged. The future climate will continue to deliver ES; in some cases, ES are increased and, in some cases, decreased-mostly decreased compared to historical supply. Along with climate change, another important factor is LULC transformation, which alters the provision of ES. Thus, mapping and evaluation of ES are crucial for future land use plans and for strengthening the capability of various services and facilities, directly or indirectly associated with ES.
Provisioning services are often readily appreciated by the public as they have direct market values. However, other ES (regulating services, supporting services, and cultural services) are often not prioritized in decision making as these services do not hold immediate monetary value. To emphasize the significance of regulating services in the production and sustainability of other ES, we focused on four major regulating ES, namely, WY, soil loss, nitrogen export, and carbon storage. The provision of these services was compared on sub-basin scales and compared based on LULC data of two time periods, namely, 2000 and 2010. We also computed the ES provision on future climate scenarios with the baseline LULC data of 2010 and downscaled climate data for 2010-2099.
The LULC change analysis between the periods 2000 and 2010 shows that the conversion to built-up and agriculture areas from other LULC is significantly high on the Bagmati Basin. This is because of increased anthropogenic activities for food sustainability and economic upliftment. These conversions come at the huge cost of increased WY, soil loss, loss of carbon sink, and increased nitrogen export to water resources. Meantime, the concept of a Payments for Ecosystem Services (PES) scheme is being promoted in the basin with the realization of the need for ecosystem conservation [55]. This scheme has increased protected areas and community forests in the basin and has a twin objective of promoting ecosystem conservation and development earnings [55]. Likewise, climate projection shows increased precipitation and temperatures for the future period. The mean annual average precipitation in the case of 2010 was lower than that of 2000 and as WY and soil loss are highly affected by precipitation amount, its magnitude was decreased in 2010 compared to 2000. Though the overall WY was decreased in 2010 compared to 2000, the WY on sub-basin 2 increased. Sub-basin 2 faced a higher increment of built-up areas and, consequently, evaporative loss and water retention were reduced, causing an overall increase of WY. This phenomenon was also observed by another study [56] when computing the impacts of LULC on WY provision. The prediction of WY on future climate scenarios shows an overall increase in all sub-basins in consecutive periods S1 to S3 under both RCP 4.5 and RCP 8.5 scenarios (Table 9). Although precipitation determines the amount of water provided by nature, LULC plays a significant role in determining the amount of water that flows as runoff or is retained as storage [57]. Owing to this phenomenon, urbanization has induced frequent urban floods in Nepal and the lack of sufficient water retention structures has increased flood peaks and flood volumes.
In terms of ES, changes in WY caused by LULC changes have two major effects, namely, contribution to the water available for consumption and/or increased flood risks during storms [58]. Infrastructural capability for water collection, delivery, and treatment may work temporarily to counter increased WY. In the long term, however, management should consider human-nature interactions to avoid unintended environmental and socio-economic consequences which are caused as a result of rapid and large-scale development [58]. Several studies have recommended ecosystem-based "green" infrastructures such as wetlands, well-structured soils, and forest patches to enhance water storage and flood regulation [59]. The mapping of WY and its tendencies in future climate scenarios can hence provide an outline for such sustainable land use plans to mitigate flood risks and water scarcity problems.
Although soil loss due to erosion and sedimentation are natural processes in a healthy ecosystem, excess loss mainly due to changes in LULC is a threat to the security of water, food, and energy [60].
Soil loss is observed to be highest from agricultural land use in the study basin. Agriculture is the major activity sustaining the economy of Nepal. Thus, the traditional agricultural system demands improved methods to sustain topsoil-containing organic matter and essential plant nutrients that are otherwise removed from the soil during erosion. Likewise, the predicted increment of soil loss in the future climate scenarios also indicates the necessity for improved farm management practices that incorporate the management of increased WY. In addition to the stress on food production, soil erosion also reduces water and nutrient retention, biodiversity, and water resources/downstream power generation [61]. Therefore, for sustainable management of soils and related services, spatial mapping of soil erosion on different land use and future scenarios can provide guidelines and mark areas for immediate actions to mitigate the loss.
The loss of carbon storage is another significant phenomenon observed in the study basin with a reduction in forest areas and increases in agricultural and built-up areas. Terrestrial-based carbon storage and sequestration are directly affected by LULC. The concentrated anthropogenic activities on the upper part of the basin, especially in sub-basin 2, contribute to the lowest carbon storage and highest nitrogen export per ha of land in this sub-basin. Most developing nations have serious problems with the degradation of forests and soils which have crucial implications for changing the global C pools [62]. The information on carbon storage and fluctuations can help land managers to choose among sites for protection, harvest, or development. Furthermore, these maps can support multiple decisions by governments, NGOs, and stakeholders to offer incentives to landowners in exchange for forest conservation. Sub-classification of land use and consideration of altitudinal variation in various carbon pools produces more accurate results. Poudel et al. [63] have found higher carbon stocks at the higher altitudinal gradient in the study conducted in the Panchase Conservation Area in Nepal. The study indicated that due to human disturbance, carbon stocks were low at low altitude. Temperature and precipitation also have significant effects on the carbon pool of various biomass [64]. Thus, a detailed study with climate and elevation variation helps to produce an accurate estimate of carbon value.
Similarly, with increased agricultural activities, and reduced vegetation covers, nutrient retention diminishes and, hence, the amount of nitrogen entering the river network/water resources increases. This severely increases water contamination, thereby affecting human and aquatic health. Additionally, because of increased WY in future climate scenarios, nitrogen export is expected to increase in all sub-basins. In this scenario, spatial information on nutrient export and areas with the highest filtration can help land use planners to integrate the contribution of ecosystems in order to mitigate water pollution.
To understand the temporal change of ecosystem services in 2000 and 2010, the minimum value was set to 0 and the maximum value was set to 1 for each ecosystem service in two periods. This shows a comparative provision of ES in all sub-basins in the 2000 and 2010 periods which depicted that sub-basin 2 had the lowest combined ES provision in both 2000 and 2010 ( Figure 12). Thus, this indicates the need for urgent policies in order to restore the services and to promote its sustainability. Likewise, the computation of ES with the projection of climate data in the 2010-2099 period shows increased WY, soil loss, and nitrogen export from the study area. The case of sub-basin 2, especially of the Kathmandu Valley area, can be referenced for the planning of other emerging cities. As Nepal is prioritizing decentralization and focusing on the development of many other smart cities on the outskirts of the Kathmandu Valley and other parts of the country, studies like this can present risk analyses and help decision making that prioritizes the conservation of ES for the optimum utilization and preservation of natural capital. of the Kathmandu Valley area, can be referenced for the planning of other emerging cities. As Nepal is prioritizing decentralization and focusing on the development of many other smart cities on the outskirts of the Kathmandu Valley and other parts of the country, studies like this can present risk analyses and help decision making that prioritizes the conservation of ES for the optimum utilization and preservation of natural capital.

Limitations of the Study
The InVEST tools used for the evaluation of ES have their own modeling limitations. The InVEST WY Model is based on annual averages and neglects the extremes and the spatial distribution of LULC. The InVEST Carbon Storage Model assumes that each LULC is at the fixed carbon storage level and the fluctuation on carbon storage is only due to a change in LULC from one type to another. Due to data unavailability in the study area, tabular values were obtained from the InVEST manual. The study used reference data for various LULC types from InVEST user guidelines and the variation in the amount of carbon storage in different carbon pools due to elevation have not been acknowledged in the study at this stage. In addition, for climate data, uncertainties remain from climate change data itself and from the selected methods of downscaling and bias correction. Nevertheless, careful attention was given to the preparation of data and modeling. Despite the model and data limitations, the results from the study provide an overview of a general tendency of the provision of ES, including fluctuations with a change in climate and LULC. It is, thus, expected to help in decision making and scenario analysis. The results of this study could be improved if ground observation data are available for the accurate analysis of InVEST models.

Conclusions
This study attempted to evaluate ES and their fluctuations with LULC and climate change on the Bagmati Basin of Nepal. The study first assessed ES based on the 2000 and 2010 LULC map and then with the baseline LULC map of 2010 and projected climate data from the MME of 12 GCMs, ES provision on the future period was estimated. The overall provision of combined ES on sub-basin 2 was lowest in 2000, 2010, and on all future climate scenarios. In addition, the provision of ES was

Limitations of the Study
The InVEST tools used for the evaluation of ES have their own modeling limitations. The InVEST WY Model is based on annual averages and neglects the extremes and the spatial distribution of LULC. The InVEST Carbon Storage Model assumes that each LULC is at the fixed carbon storage level and the fluctuation on carbon storage is only due to a change in LULC from one type to another. Due to data unavailability in the study area, tabular values were obtained from the InVEST manual. The study used reference data for various LULC types from InVEST user guidelines and the variation in the amount of carbon storage in different carbon pools due to elevation have not been acknowledged in the study at this stage. In addition, for climate data, uncertainties remain from climate change data itself and from the selected methods of downscaling and bias correction. Nevertheless, careful attention was given to the preparation of data and modeling. Despite the model and data limitations, the results from the study provide an overview of a general tendency of the provision of ES, including fluctuations with a change in climate and LULC. It is, thus, expected to help in decision making and scenario analysis. The results of this study could be improved if ground observation data are available for the accurate analysis of InVEST models.

Conclusions
This study attempted to evaluate ES and their fluctuations with LULC and climate change on the Bagmati Basin of Nepal. The study first assessed ES based on the 2000 and 2010 LULC map and then with the baseline LULC map of 2010 and projected climate data from the MME of 12 GCMs, ES provision on the future period was estimated. The overall provision of combined ES on sub-basin 2 was lowest in 2000, 2010, and on all future climate scenarios. In addition, the provision of ES was observed to be decreasing in all other sub-basins. Outcomes like increased WY, reduced carbon storage, increased nitrogen export, and soil loss suggest that immediate actions are required from policy makers for the sustainable management of natural capital. The ES are interrelated and the absence of adequate designs to sustain one can hamper other provisions as well. The availability of spatially explicit information on ecosystems and their interrelated services serves for the prioritization of ES into policy and decision making. The maps produced as an outcome of this study can help land use planners, government organizations, and concerned stakeholders to recognize areas where the ecosystems are produced and help in the decision making for low impact development maintaining ecological balance and economic goals. Further studies on national scale future LULC scenarios could estimate accurate ES provisions integrating national priorities.
Supplementary Materials: The following are available online at http://www.mdpi.com/2071-1050/12/19/8186/s1, Table S1: list of stations with rainfall and temperature data along with their location and data period. The data has been duly collected from Department of Hydrology and Meteorology (DHM) Nepal.