Quantifying water-related ecosystem services potential of the Kangchenjunga Landscape in the eastern Himalaya: a modeling approach

Ecosystem services (ES) provide direct and indirect benefits to human well-being as a result of ecosystem structures and functions (de Groot et al. 2010; Burkhard et al. 2012). These ES can be grouped as provisioning, regulating, cultural, or supporting services (MA 2005). Burkhard et al. (2014) highlighted that the ecosystem service potential of an area differs from the real ecosystem service gained by the society. The ES potential refers to the maximum capacity to supply selected ES, based on their functioning and the demand for these potential services from society (beneficiary) converts them into the real ES or ecosystem flow. In other words, ES potential is the service provided by the ecosystem before it is used for benefit. The ES potential depends on the biophysical data whereas the ecosystem flow requires information on socio-economy. It is therefore important to distinguish between the potential and the actual flow of the ES.

In the past, the focus was on valuing ES through the economic lens (Chan et al. 2006; Naidoo & Ricketts 2006; Anderson et al. 2009); however, there is a shift in research, to understand the process and their continuous supply of benefits (Francesconi et al. 2016) for their sustainable use. Recently, the concept of ES has gained an enormous attraction from the decision-makers and environmentalists, to integrate them into development plans and policies (Müller & Burkhard 2012; Shrestha et al. 2021).

Water plays a significant role not only by providing direct benefit for consumption but even a key role in maintaining and regulating various other ES, such as energy and food production, habitat for fish, water cycle, and cultural and recreational activities. Therefore, water is an essential element of the ES concept and it is defined as the ES provided on a larger scale, for example, at a watershed level (Schmalz et al. 2015). Quantifying ES is a complex process. The method for measuring ES value depends upon the objectives, research questions, data availability, and domain of the ecosystem. For example, a provisioning service of a forest area to provide timbers could be evaluated in terms of monetary value. However, cultural services, which can only be measured qualitatively, are highly subjective and can have only a local or regional value. Nevertheless, the development of either quantitative or qualitative methods to measure ES is recognized as a crucial step toward its management (de Groot et al. 2010; Boyanova et al. 2014).

A hydrological model can assess various water balance components that can be used as an indicator for various WES (Nedkov et al. 2015). These indicators can serve as a measure of environmental phenomena, used to represent baseline environment conditions or changes (Müller & Burkhard 2012). The selection of freshwater-related indicators from the hydrological models is crucial and it should have a relevant relationship to the ES (Kandziora et al. 2013). Mapping of various WES comes in handy to identify linkages across many realms of water resources management (Vigerstol & Aukema 2011; Liu et al. 2013). Maps can help visualize data spatially and temporarily to understand patterns and show potential resource limits which can support decision making. This information is also vital for understanding future changes in the ES under climate change scenarios. The selection of the model(s) depends on the input data availability, its resolution, climatic zones, location of the study area, and the output the model can process. For example, the J2000 model will be appropriate for estimating snowfall and snowmelt in a monsoon-dominated Himalayan region rather than the InVEST model.

Watershed management and biodiversity conservation often exceed geopolitical boundaries (López-Hoffman et al. 2010; Gurung et al. 2019; Kotru et al. 2020). The landscape concept has gained worldwide attention (Hamilton & McMillan 2004; López-Hoffman et al. 2010; Sayer et al. 2013; Dallimer & Strange 2015; Reed et al. 2016; Gurung et al. 2019; ICIMDO 2019) as being vital for both long-term resource management and equitable development across the border (Kotru et al. 2020). One such landscape with a history dating from 1997 is the Kangchenjunga Landscape (KL). In 2002, ICIMOD initiated the Kangchenjunga Landscape Conservation and Development Initiative (KLCDI) with stakeholders from Bhutan, India, and Nepal. The landscape is located in the eastern part of the Hindu-Kush Himalayan (HKH) region (Xu et al. 2019). Covering approximately 4.3 million km² area, the HKH region is characterized by rich ES and diverse biodiversity. Also known as the ‘water tower of Asia’ it is estimated that the total renewable annual water availability of the HKH region (eight countries) is 7,745.5 km³ (FAO 2016). Of the total renewable water, 20.62% or 1,597.8 km³ is used annually for agriculture, industrial, domestic, hydropower generation, and other purposes (Scott et al. 2019). The ES from the HKH region provide a livelihood to more than 240 million people living within the HKH region and an additional 1.9 billion living downstream (Sharma et al. 2019).

This study aims to estimate the water-related ecosystem services (WES) potential in the KL, using a hydrological model. The rivers originate within the KL and flow towards the south, where it meets the larger Ganges–Brahmaputra River system. The landscape is home to more than seven million people, 4,500 species of plants, 169 mammals, and 618 species of birds (Gurung et al. 2019). Rich in water resources, the KL supports larger ecosystem services including hydropower, agriculture, and tourism and generates a significant amount of revenue for the KL member countries. A large number of tourists visit the landscape to see hill stations, tea gardens, monasteries, lakes, rich wildlife, and the third highest mountain, the Kanchenjunga. This is the first study that we know of to utilize the process-based cryospheric hydrological J2000 model, to quantify the WES potential at a landscape level. From the model outcome, we looked at the spatial and temporal water balance components to quantify various WES.

The KL stretches along the southern side of the third highest peak, Mount Kangchenjunga (8,586 masl), and covers an area of 25,081 km². The KL spreads across parts of eastern Nepal (21%), Sikkim and West Bengal of India (56%), and the western and south-western parts of Bhutan (23%). Mount Kangchenjunga, which lies at the heart of the landscape, is considered sacred by the Tibetan, Sikkimese, and Kiranti communities. The elevation of the landscape ranges from 40 masl in the south to above 8,000 masl to the north (Figure 1). Five major rivers, Tamor (Nepal), Teestaand Jaldhaka (India) and Torsa, and Wang Chu (Bhutan), flow through the KL and meet the mighty Ganges and Brahmaputra river system.

The landcover area of the KL was analyzed based on the ESA map of 2010. The largest part up to 36% of the landscape is covered by forest areas (Evergreen and Deciduous). The cropland area covers 24% of the total area, and it is predominantly located in the southern part of the KL. Grassland and shrubland each cover an area of 8%. Roughly 6% of the area is permanently covered by snow and 3.5% by glaciers in the high altitude areas from where the Tamor and Teesta rivers originate. The soil type of the KL could be roughly divided into five types based on the Harmonized World Soil Database (Figure 1). Based on the USDA texture classification, the landscape has a soil classification of loam. Loam soil is mostly composed of sand, silt, and a small amount of clay.

Figure 2 provides a snapshot of the methodological approach used for this study. The J2000 hydrological model was used to simulate the various eco-hydrological variables to estimate WES. The J2000 model is a process-oriented distributed model based on hydrological respond units (HRUs). The HRU is the smallest spatial unit of the model made up of similar land uses, soils, and slopes within a subbasin based upon user-defined thresholds (Flügel 1995). The model can also replicate cryospheric processes such as snow and glacier melt.

The KL model (KL-J2 K) covers the area beyond the boundary of the KL including the upstream of Torsa and Wang Chu River Basins, and the downstream of Tamor basin up to the Majhitar hydrological station for calibration (Figure 3). The model KL-J2 K is composed of 21,713 HRUs, of which 18,133 HRUs lie within the boundary of the KL. This brings the average area of the HRU to be 1.43 km².

Only two rainfall stations and four temperature stations on the Nepal side, and five rainfall and five temperature stations on the Bhutan side in the KL were available during the period of study (Figure 3). Two hydrological stations, namely Majhitar (Tamor) and Chukha (Wang Chu) lie in the KL, which were used for model calibration and validation. Since the in situ (station) data were only available on the Nepal and Bhutan side of the landscape, the APHRODITE rainfall data were used to compensate for precipitation data scarcity.

The model was set up and calibrated in the following three steps:

• Model parameter values transferred from a published Tamor J2000 model.

• Calibration of the APHRODITE rainfall data to simulated discharge.

• Final calibration using snow cover and discharge data.

Since the Tamor basin is also part of the KL, model parameter values from a published calibrated Tamor J2000 model (Tamor-J2 K) were transferred to the KL-J2 K (Nepal et al. 2017b). Transferring model parameters to a nearby basin with similar biophysical characteristics has been a common practice in places where there is a lack of data (Nepal et al. 2017b; Shrestha et al. 2017a). Nepal et al. 2017b have conducted a detailed study to show model transfer suitability using the J2000 models. Therefore, this study uses the parameters from the published Tamor-J2 K to perform further analysis.

Finally, the model was fine-tuned to simulate the snow cover area of the entire landscape and discharge at the two aforementioned hydrological stations. ICIMOD's MOYDGL06 data set, developed by combining images from Terra and Aqua MODIS 8-day snow cover, was used to calibrate the model (Muhammad & Thapa 2020). Widely used goodness-of-fit statistics like r², Nash–Sutcliffe efficiency (NSE), Pbias, and root mean square error (RMSE) were used to evaluate the performance of the model (Shrestha et al. 2017b, 2018).

The indicators of ES potential for the KL were selected based on the output available at the HRU and a river basin level. Table 1 provides a list of ES potential indicators used in the study, for analyzing provisioning and regulating services (Burkhard et al. 2014; Schmalz et al. 2016) along with the corresponding output variable from J2000. Provisioning services are the products that are obtained from nature for direct consumption or use. Here, the water yield at the basin outlet for each of the five rivers and a combined water yield at a landscape level are analyzed. Water yield refers to water available (m³/s) at the outlet of each river segment. The information is analyzed at seasonal and annual resolutions. This information is quite handy for analyzing water availability for domestic uses, irrigational, hydropower generation, industrial, environmental flow, and even recreational and tourism services (Francesconi et al. 2016). The contribution of surface and groundwater (interflow and baseflow) flows is also estimated.

The amount of precipitation, snowmelt, evapotranspiration (ET), soil moisture (SM), and infiltration are important indicators of regulating ES. This information is available at the HRU level. Annual and seasonal precipitation, snowmelt, and infiltration are important regulating service indicators that govern the amount of water available in a landscape. ET and SM are other important indicators that have a direct impact on agro-ecosystem production. Higher ET and low SM can have an adverse impact on vegetation, and they can impact on habitat and ecosystem of the landscape. These ES potentials are presented spatially by aggregating annual and seasonal average values over the period of 2000–2015.

To check the accuracy of the APHRODITE data, they were compared against the few available precipitation stations. Figure 4 shows the average monthly comparison of precipitation from the APHRODITE and measured rainfall stations. We used the APHRODITE grid, which overlaps with the near station for the comparison. The analysis suggests that the APHRODITE data show a lower amount of rainfall during the monsoon season. Other studies also observed lower rainfall values in APHRODITE than in observed stations. For example, 13–25% less rainfall was found in the Tibetan Plateau by Tong et al. 2014 and, up to 200% in the upper Indus basin was found by Immerzeel et al. 2015. Similarly Ji et al. (2020) found that when comparing against the in situ data at the Yarluzangbu–Brahmaputra River basin, the APHRODITE data had a relative bias of up to 40% in some stations.

To compensate for the underestimation in the APHRODITE precipitation, the value was increased gradually from 10 to 60% to simulate the discharge at Majhitar and Chukha stations. The efficiency criteria (NSE and BPAIS) were tracked for each simulation, which helped to find out the best combination to satisfy the water balance at the daily level. It was found that by increasing the APHRODITE value by +40%, the model performance increases to a maximum in terms of NSE and least bias at both Majhitar and Chukha stations (Figure 5).

After the rainfall adjustment, the KL-J2 K was further calibrated manually first against the snow cover area, and then to discharge data at Majhitar (Tamor) and Chukha (Wang Chu) stations. The initial model parameters were transferred from the Tamor-J2 K as discussed earlier. Figure 6 compares the observed and simulated annual snow cover area for the period 2002–2010. The annual snow cover area was captured well by the model except for the years 2006 and 2010. The model was able to simulate the daily snow cover area of the entire KL, with an accuracy of r²=0.47, and RMSE of 1,452 km², and the annual snow cover area of r²=0.88, and RMSE of 402.1 km². An earlier attempt at snow cover validation in the western part of the KL, i.e. Tamor basin, with the GR4JSG model was also 0.46 (Nepal et al. 2017a).

The total simulated period for discharge was divided from 2000 to 2009 (10 years) as a calibration period and 2010–2012 (3 years) for Tamor and 2010–2014 (5 years) for Chukha as a validation period, with 1999 as a warmup period. The model performs well during both periods at both stations (Figure 7). During the calibration period, the model performed with an NSE of 0.83, r²=0.85, and Pbias −11.96% at Majhitar and with an NSE of 0.83, r²=0.83, and Pbias 4.5% at Chukha stations. Similarly, during the validation period, the NSE, r², and Pbias were 0.86, 0.88, −17.55%, at Majhitar and 0.75, 0.81, 18.67%, at the Chukha station.

Four distinct seasons could be seen in the KL, namely winter (DJF), spring (MAM), summer (JJAS), and autumn (ON). The total annual precipitation and the mean average temperature of the landscape are estimated at 2,790 mm and 12 °C, respectively. Nearly 90% of the precipitation falls during the months from April to September. The month of July experiences the highest precipitation with 650 mm, whereas December has the least precipitation of less than 10 mm. During the winter seasons, the KL has collective precipitation of 60 mm and the minimum temperature could reach up to −2 °C. During the summer seasons, the landscape-wide average maximum temperature can reach up to 18 °C. Precipitation in the KL is heterogeneously distributed both temporally and spatially (Figure 8). The precipitation amount increases towards the southern flats area. Of the total precipitation that falls in the KL, 4% (112 mm) falls as snow annually. The monthly snowfall percentage however varies significantly. During the winter season, snowfall percentage ranges from 20 to 30%.

Figure 9 shows the spatial map of the actual ET distribution for each month. The lowest landscape-wide monthly ET of 46 mm is observed during December, whereas the month with the highest ET at 119 mm is observed during May. With an elevation range from 40 to more than 8,500 masl, the ActET also varies considerably. The ET can reach more than 20 mm/day during the hottest month in the southern part, whereas it can reach up to nearly 0.1 mm/day during the coldest months in the higher elevation. The average annual actual ET of the landscape is estimated at around 995 mm. Figure 10 provides information on annual spatial SM in the KL. SM is the water stored in the soil and governs the exchange of water between land and atmosphere via ET, and hence is an essential element for crop production. SM is affected by precipitation, temperature, and type of soil. The SM of the landscape varies from 330 to 0.15 mm/day. The SM is found higher in the mid-hill, in areas covered by forest.

Roughly 70% of the total precipitation is converted to discharge annually (Figure 11). Separating the runoff into three components reveals that the surface flow covers 65% of the annual runoff, whereas interflow (flow within the soil) and the baseflow (from deep aquifers) cover 21 and 14% of the total runoff, respectively. Surface flow, however, covers a significant amount of total runoff that is only concentrated during the monsoon period. The interflow and baseflow are relatively constant annually, and they provide regular water in the rivers during the dry season. The glacier and snow cover area in the KL, covering less than 10% of the total area, are not only important for climate regulation but also supply annually 11% of the total water.

Figure 12 provides a comparison of regulating ES potential of the five basins in the KL. Jaldhaka basin, although covering the smallest area (Table 2), receives the highest amount of precipitation. 100% of the total precipitation in the basin falls in the form of rain. On average, the basin experiences 4,600 mm of rainfall in a year. The actET, which correlates with the temperature of an area, is highest in Jaldhaka than in the other basins. The actET in the basin is 1,635 mm which is 36% of the total precipitation. Teesta, Torsa, and Wang Chu, the three largest basins of the KL, are spread from the north to south (Figure 1), and as a result a large spread in the box plot of precipitation and actET is seen for the basins (Figure 12). On average Torsa, Teesta, and Wang Chu experience 2,900, 2,500, and 2,000 mm of precipitation in a year, respectively. The smaller Tamor basin experiences higher rainfall of 2,070 mm compared to the larger Wang Chu. Tamor has the highest percentage of snowfall, which is 13% of the total precipitation. The snowfall percentage for Teesta, Torsa and Wang Chu lies in the range of 3–6%. Tamro has the highest ice melt of 200 mm, followed by Teesta with 126 mm per year. For Torsa and Wang Chu, the ice melt values are 63 and 60 mm, respectively. This shows that the snow and glacier melt has a significant role in generating discharge in the rivers. Almost 40% of the total precipitation in the Wang Chu basin is lost in the form of ET. The annual SM of all five basins is similar with the value ranging from 115 to 130 mm per year. The SM in all the basins is relatively similar as shown by a narrow spread in the box plot.

Table 2 provides the annual average and seasonal water yield for each of the five rivers, and at the landscape level. The total average annual water yield of the KL is estimated to be 1,908 m³/s (112.7 million m³). The annual discharge from the Tessta River contributes 28% of the total yield in the KL followed by Torsa and Wang Chu with 16 and 14%, respectively. With an average annual discharge of 196.6 m³/s, Tamor contributes 10%, whereas Jaldhaka contributes 11% of the total yield. Other smaller rivers contribute roughly 22% of the total yield. The contribution from snowfall and ice melt in the Jaldhaka basin is zero (Figure 12), implying that 100% of the water availability is generated through rainfall. The temporal distribution of water yield at both landscape and basin-wide are heterogeneously distributed. Water yield during summer, which is also the monsoon period in the region, can reach up to an average of 4,416.8 m³/s. The water yield in the landscape decrease to 1,136 m³/s during the autumn or the post-monsoon periods, and the lowest water yield is seen during the winter season with an average discharge of 360 m³/s.

The results from any model are always associated with a certain degree of uncertainty. For a hydrological model, the uncertainty is usually related to either the input data, the selection of the model and algorithms, or the combination of both. In the case of the KL, precipitation data were a challenge. This challenge was mostly bypassed using the APHRODITE data, which has been proven suitable in the Himalayan region by many studies (Ji et al. 2020; Thapa et al. 2020). As with any model, it is especially important to use as accurate and representative information as possible, to gain meaningful results. We have compared the APHRODITE data against the available observed rainfall stations in Nepal and Bhutan. We found that the APHRODITE data are consistently underestimating in most of the stations during the monsoon period. The average difference between the observed and the APHRODITE data was around 37% (56% in Nepal and 18% in Bhutan). Therefore, a landscape-wide multiplicative correction factor was used to correct the APHRODITE data before calibrating the model. We performed the experiment to increase the APHRODITE data from 10 to 60% and tried to simulate the discharge. Other studies have also shown such underestimation, especially at higher altitude.

A well-established, J2000 hydrological model was used to represent the water balance of the KL. The model is suitable for calculating not only snowmelt but also glacier melt contribution (JAMS Models 2020). Hence, it is better suited to represent the basins of the Himalayan region. Spatial information of 90 m resolution was used to set up the model, which produced 21,713 individual HRUs. This size of HRUs is suitable to simulate the ecosystem service information at a landscape level.

Although the model performed well in simulating the discharge data at Majhitar and Chukha stations for Tamor and Wangchu rivers (Figure 7), the performance at the other three rivers could not be evaluated due to data availability. To make the model robust, it was also calibrated against the satellite-based snow cover information at a landscape level. In addition, since snow is an integral part of the landscape, it is also equally important that the model should be able to simulate snow components. By calibrating the snow cover area and streamflow, we believe that snowfall and snow/ice melt are fine-tuned to some extent and thereby help us to estimate the contribution of snow to the total water yield in the KL. The model was able to simulate the daily and annual snow cover area with the accuracy of r²=0.47, RMSE of 1,452 km² and r²=0.88, RMSE of 402.1 km², respectively.

Annual result is less relevant in answering water-related ES in areas where the hydrology events are heterogeneously distributed. It would therefore be more representative to produce the information at the monthly and seasonal level. Around 90% of the precipitation falls during the monsoon periods, from April to September, with the southern region receiving the majority of the precipitation. Precipitation in the form of snow is an important regulating service in the KL. They fall in the northern higher elevations. Although it accounts for only 4% of total precipitation annually, it is much more prominent during the autumn and winter seasons. Snowfall during these seasons stores water, that slowly melts during the thaw seasons, providing valuable water resources, especially in the northern area. The annual average ET of the landscape is estimated to be around 995 mm which is around 35% of the total precipitation. ET is also termed green water, which is the water used by vegetation for its growth. This green water is a good strong indicator of ecosystem productivity. In higher elevations of the KL, where rain-fed agriculture is dominant, vegetation solely relies on green water.

SM is a function of soil property, vegetation cover, and atmospheric condition and it is one of the major requirements for sustaining biomes. Although in the J2000 model there are two soil storage components: Middle Pore Storage (MPS) and Large Pore Storage (LPS), for this study, water stored in the MPS has only been considered as SM. The moisture in the MPS represents the field capacity, whereas LPS represents water contribution from the shallow aquifer. The input for the SM is from the rainfall and snowmelt, whereas ET is responsible for soil dryness. Water in the MPS is used by the plants for the ET process, hence providing a vital regulating service for ecosystem productivity. Around 60% of the total area in the KL is covered by mixed forest and cropland. The trees and crops that are sustained because of SM provide habitat for animals and timber, fuel, and food for humans. Wheat, rice, maize, and tea are usually grown in the region. Another important cash crop in the KL is the large cardamom. Large cardamom is one of the primary export commodities in Bhutan (18.6 m USD), India (11.6 m USD), and Nepal (34.5 m USD) (ICIMOD 2019). These crops require shade trees to receive both light and shade, and they require ample moisture to survive. Therefore, the amount of SM directly affects water availability, crop yield, economy, and ecosystem productivity.

About 70% of the total precipitation in the KL is converted to blue water. The total average annual water yield of the KL is estimated to be 112.7 million m³ (1,908 m³/s). This figure is estimated by summing the yield of individual river outlets. Also termed as blue water, it represents both surface and groundwater that can be used for irrigation, power generation, industrial, domestic, and recreational uses. Teesta Basin alone has about 47 reservoirs in various stages of development (Rahaman & Mamun 2020) used for hydropower, and irrigation. Similarly, 14 HPP with the capacity to produce 109 MW of energy is in operation (DOED 2021) in the western part of the KL within the Nepal border. Plans are in discussion to divert water from Tamro basin to the Morang district for irrigating 114,000 ha of land in addition to generating 732 MW of energy (DWRI 2019). Around 30 species of fish consider the Tamron River as their home with most of them having economic value as food, medicine, and/or aesthetic value (Shrestha et al. 2010). Hydropower is the highest revenue generator, contributing to 80% of the export earnings in Bhutan (Alam et al. 2017). Wangchu River basin in Bhutan plays an important role in its economic and social welfare. It supports two of the largest hydropower plants; Chukha (336 MW) and the Tala (1020 MW) (MOES 2018). Chukha hydropower plant alone contributed over 30% to Bhutan's total revenue in 2006.

The water yield can be further divided based on contribution from ice melt (from snow and glacier), and direct runoff, or based on surface and sub-surface runoff. Snow and glacier melt provides flow regulating services during the dry season. They account for 14% of the total water yield in the KL, despite the total area coved by snow and glacier being less than 10%. Four out of the five major rivers originate and rely on the melting of the snowpacks and glaciers. These ice melts are important to ensure that the rivers remain perennial, and provide water supply services throughout the year. These rivers flow downstream and are an important tributary of the Ganges–Brahmaputra River basin.

Groundwater also plays a significant role in providing freshwater throughout the landscape. In this study, the interflow and the baseflow components are used as an indicator of groundwater. The interflow and the baseflow are the contribution from the shallow and deep aquifer to the streamflow, respectively. The combined interflow and baseflow account for nearly 35% of the total water yield annually. Although these flows are constant throughout the year, they are much more essential during the dry seasons. Results suggest that during the dry seasons, the contribution from interflow and baseflow is nearly 100% of the total water yield in the rivers. These components of flows are a function of LULC and soil type. In general, areas covered by vegetation are suitable for infiltration and act as a recharge area. Conversion of such recharge areas to a non-permeable area, like urban build-up, can have an adverse impact on seasonal water availability in the landscape.

This study attempts to quantify the WES in the KL using suitable indicators. The data produced provide relevant information regarding the hydrological components and their importance in sustaining the ecosystem productivity of the landscape. The information is produced in a spatio-temporal pattern, which can be used by the stakeholder for further discussion and planning of sustainable land management.

This study is a part of the Kangchenjunga Landscape Conservation and Development Initiative of ICIMOD. This study was supported by the core funds of ICIMOD, contributed by the governments of Afghanistan, Australia, Austria, Bangladesh, Bhutan, China, India, Myanmar, Nepal, Norway, Pakistan, Switzerland, and the United Kingdom. The authors would also like to acknowledge the Department of Hydrology and Meteorology, Nepal, and the National center for hydrology and meteorology, Bhutan for providing data for the study. We also acknowledge the help of Dr Tashi Dorji for helping us acquire data from Bhutan. Finally, the authors thank Kabi Raj Khatiwoda, Dr Nakul Chettri, Nisha Wagle, Pradeep Dangol, Dr Randall Ritzema, Sanita Dhaubanjar and Saurav Pradhananga for their technical support and valuable suggestion.

The views and interpretations in this publication are those of the authors and are not necessarily attributable to ICIMOD.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.