Simulating glacier mass balance in the cross-border Poiqu/Bhotekoshi Basin, China and Nepal

In many mid-latitude mountain regions, glaciers act as a temporal water storage such that new snow falls to nourish the glaciers in winter and then melts in summer (Jansson et al. 2003). Normally, the glaciers in these regions can keep a balance between the melting and nourishing, however such balance is at risk due to the climate change. Recently, glacier shrinkage has become increasingly serious and has drawn much attention from the public. The most crucial concern is that the shrinking glaciers are likely to cause water loss and stream-flow variability, particularly during the dry seasons of the year (Braun et al. 2000; Collins 2008; Moore et al. 2009). At transnational level, the concerns regarding the mid-latitude glacier shrinkage mainly concentrate on the accessibility of water resources, which has caused disputes and even conflicts among countries. Poiqu/Bhotekoshi basin, a mid-latitude transnational region between China and Nepal, is currently confronting the imbalance between glacier recharging and discharging (Fujita & Ageta 2000). A quantitative assessment of the glacier loss is urgently needed in order to support environmental decision making.

In this study, a hydrological conceptual framework is considered to support the quantification of the glacier loss. The conceptual model used by this study is different from the physical model. The physical model is based on the hydrological process research using experimental methods. However, some measured data of cross-border basins are very difficult to obtain. The conceptual model generalizes the physical basis of the basin, and then combines the hydrological empirical formula to approximate the flow process of the basin. The basic idea is to model the runoff during the melting season as compared to the rainfall or snow nourishment such that the GMB can be estimated. The main challenges of this study come from several aspects. First, the existing hydrologic models that are mostly developed for estimating the contribution of melted ice sourced from the mountain snow peak, which could not provide a straightforward inference for estimating the runoff caused by mid-latitude glacial melting. It is worth noting that many scholars have made great progress in the study of glacial runoff at very high altitudes in recent years (Zhu et al. 2016; Yin et al. 2017), especially in the estimation of runoff using the improved HBV model at a glacierized alpine catchment (Wang et al. 2019). Second, the runoff in high elevation upstream basins can be contributed by either ice or snow, and therefore an analytical separation is required to identify the flows from the melting glaciers. In addition, the bad transportation conditions in the study area led to the rarity of in-situ hydrological data. Despite the dataset being adequate in the context of temporal information, only 12 weather stations were sparsely distributed over the basin. To resolve the abovementioned concerns, a runoff model that provides separate estimations for ice and snow and is capable of handling ‘small’ datasets is significantly necessary (Kulkarni et al. 2002; Zhao et al. 2013; Sun et al. 2015).

A large number of previous studies have shown that the methodological development in glacial hydrology has gone through three stages, namely empirical statistical modeling, analytical modeling, and numerical modeling (Chen et al. 2003). The first two stages appeared in the 1970s and have no obvious temporal order in the application to estimate glacier runoff. Since the early 1990s, simple statistical models started fading in the applications of simulating glacier runoff. On the other hand, the conceptual glacier hydrological model was in development on the basis of formulating the physical process of the glacier runoff. Until the late 1990s, Arnold et al. (1998) applied the numerical modeling technique ‘distributed hydrological modeling’ to study the glacier hydrology. Distributed hydrological modeling could truly describe and scientifically reveal the natural spatio-temporal dynamics in the context of the hydrological circle, and opened a new era of glacier hydrological modeling (Arnold et al. 1998).

Since the HBV model was successfully applied to predicting the possible change of the runoff in an inland mountainous river basin under the conditions of climate change, many variations of HBV have been developed (Bergström 1976; Kang et al. 1999; Mclntyre et al. 2013). The advantages of using HBV-based models mainly come from the fact that they can estimate the runoff from both snow and ice melting, and are also highly tolerant to a small number of training samples. Moreover, HBVs are capable of modeling the glacier runoff across different river basin regions (Bergström 1992; Jost et al. 2012). The studies of HBV have focused on the topics of sensitivity and uncertainty analysis targeted at the model parameters since the 1990s (Abebe et al. 2010; Chen et al. 2012; Li & Xu 2014; Nourali et al. 2016), the contrastive study between the HBV model and other hydrological models (Clark & Kavetski 2010; Kavetski & Clark 2010), the model application under different climate conditions (Braun & Renner 1992; Liden & Harlin 2000; Li et al. 2013), and integrating the HBV model with climate models through which the influence of climate change on the runoff can be investigated (Menzeil & Burger 2002; Gardelin et al. 2003; Hagg et al. 2007).

As one variation of the conceptual HBV model, the glacier melt runoff was rescaled by elevations in order to be adapted to simulate the Heihe mountain runoff in the north of China (Kang et al. 1999, 2002). Despite Kang et al. (1999, 2002) computing the runoff contributions from ice/snow and rain separately, their description of the runoff formation process was incomplete, for instance, the infiltration and melt water refreezing were not considered. This study added infiltration and melt water refreezing factors to the model, and a complete conceptual HBV model was proposed and then applied in the Poiqu/Bhotekoshi basin to simulate the runoff process and assess the changes of GMB.

The objectives of this paper are: (a) to establish a conceptual hydrology model that simulates the glacier melt runoff in the cross-border region between the Poiqu basin, China, and Bhotekoshi basin, Nepal; (b) to estimate the sensitivity of simulation factors for this model in a glacier-feed zone transitioning from a higher altitude to a lower altitude; and (c) to assess the change of GMB in Poiqu/Bhotekoshi basin and support the environmental decision making.

Poiqu/Bhotekoshi, a tributary of the Cauchy River in the Ganges upstream, is one of the main water resources of Nepal, and since it flows across three countries, it is an important and invaluable topic in transnational environmental management (Figure 1). Poiqu/Bhotekoshi river, which originates in Nyalam County, China (85°25′–E86°30′E, 27°30′–N28°35′N), is 117 km long covering an area of 2,018.41 km². In this river basin, the altitude changes rapidly from 4,554 to 616 m, which leads to a convergence of various climate, soil types, and landscapes. The precipitation in Poiqu/Bhotekoshi basin is spatially imbalanced; the northern Himalayas have 300–400 mm yearly rainfall, while the annual precipitation usually stays between 1,000 and 1,500 mm in the south (subtropical and tropical area). The annual rainfall in the middle Poiqu/Bhotekoshi basin can also reach more than 1,000 mm, however a large amount of annual evaporation is observed. The ecological degradation, population concentration and the unbalanced distribution of precipitation, which in turn lead to serious soil erosion, attribute a dry characteristic to the valley area. The basins are widely spread from the most southern region to Zhangmu-Nyalamu region (above 5,000 m), and the corresponding surface soil content changes from southern alluvial soil and mountain red soil to the yellow brown soil and cold desert soil. From south to north, the vegetation type changes from coniferous/broad-leaved to shrub/meadow. In this study area, the relative elevation difference is large, so the simple statistical model is not suitable for glacier runoff simulation. The introduction of conceptual glacier hydrological model HBV focuses on the physical process of glacier runoff formation, which is helpful for the accuracy of simulating runoff and material balance.

Due to the complexity of the conceptual HBV model, various factors (model input) are required, including daily average temperature, daily precipitation, monthly potential evaporation, temperature gradient, and annual precipitation gradient. In this study, the daily precipitation, daily temperature and monthly evaporation data were collected from 10 weather stations during the period 1988–2008. Among these stations, the upstream stations include Nielamu, Gumthang, Sundarijal and Zhbin, the midstream stations include Jiri, Dolalghat and Busti, and the downstream stations include Pandhera, Dobhan and Panchkhal. Runoff observations were acquired from three hydrometric stations from 1988 to 2008 (Figure 1). The meteorological data units involved in the model are different, but the standard unit values of snow water equivalent, evapotranspiration, soil storage and runoff are finally obtained through the calculation of nonlinear equations. To fill the dates with no meteorological observations, a combination of the precipitation gradient technique and multivariate regression was applied to reproduce the data. If there are missing daily rainfall data between 1988 and 2008, these two algorithms can fill the gap of missing daily rainfall data. The land use/land cover (LULC) map was obtained by using object-oriented classification based on the Anderson classification system (Anderson 1976). The accuracy of classification results is 84.5% compared with field validation. This data was used in the model setting, that is, land use data was used to calculate land use area at different altitudes in different ecosystems. As observed from the weather stations, the temperature drops at a rate of –0.61 °C/1,000 m as the elevation increases. However, precipitation has a positive relationship to the elevation: an average 25 mm precipitation increase corresponds to every 1,000 m elevation rise.

An HBV technique, considered as a semi-distributed conceptual basin hydrological model, was developed by the Swedish Meteorological and Hydrological Research Institute based on the characteristics of the basins in the Nordic region (Bergström 1976). Its advantage lies in the parameterization of the physical processes of runoff generation and confluence based on the principle of water and energy balance, which has a clear physical concept. HBV is a semi-distributed model in which the catchments could be separated into three elevations and seven vegetation zones. By using daily precipitation and air temperature data, such a model may be used to simulate the basin daily runoff (Figure 2) (Seibert & Vis 2012). The daily meteorological hydrological data from 1988 to 2008 were divided into two groups (the standard value of 1988–1998, and the reference value of 1999–2008). By using the Thiessen polygon method, the spatio-temporal interpolation of continuous daily precipitation value was conducted based on the precipitation data of six stations (Nyalum, Gumthang, Barhabise, Dhap, Nawalpur and Chautara). Then, the temperature data in two of the stations (Nyalum and Panchkhal) were used to estimate the daily average temperatures of the Poiqu/Bhotekoshi basin.

Equations (5)–(7) were proposed as a variation of the traditional conceptual HBV model. Our modifications mainly take into account the seasonal and annual variation characteristics of the glacier system. Newly developed equilibrium line altitudes were deduced from our observations and then employed to standardize the model calibration process for the regions with different elevations, such as accumulation zone and melting zone. The determination of the glacier equilibrium line height is very important for the calculation of the HBV model, because it can eliminate the calibration fuzziness of model parameters SFCF, CFR and FC (Table 1).

The daily observations were partitioned into three different time intervals: 7, 14, and 30 days. These three values are the set values of the model, reflecting the changes of runoff in 1 week, 2 weeks and one month. In accordance with the runoff observations in the study area, we defined the runoff volume, Q_obs, as follows: less than 0.1 as extreme dry year; 0.1–0.5 as dry year; 0.5–1 as normal year; 1–2 as wet year; and Q_obs > 2 as extreme wet year. These five thresholds are obtained by repeated calculation of the model, reflecting the real five grade indexes of the runoff anomaly change of Poiqu/Bhotekoshi Basin. Based on the above definition, the annual runoff data were classified into different categories. To prevent duplication, the data from closely adjacent stations were aggregated and allocated to the lower reaches of the basin. To avoid the high computational cost, the altitude was discretized into three levels: <1,800, 1,800–3,000, and 3,000–5,000 m. The land cover types were organized into two major ecosystems – glacier and forest (Table 2).

To calibrate the model parameters, we combined a genetic algorithm (GA) with a trial and error method, along with visual examination. The R_eff was used for the simulation of the standard value, and the reference value was used for the parameter setting. High quality data (1988–1998) was used for the parameter calibration in order to obtain a comprehensive model parameter (Table 2). The sensitivity of the model was tested for 15 parameters (Table 3), and it was found to be significantly sensitive to recession coefficient K₁, the snow adjustment factor and the soil field water-holding capacity, and also relatively sensitive to threshold temperature (TT), refreezing meltwater (CFR), a certain portion of the water equivalent of the snow pack (CWH), lower base flow recession coefficient (K₂) and max percolation to lower zone (PERC). In order to test the model efficiency, the calibrated parameters were used to conduct the non-parameter calibration. In this study, the data from 1999 to 2008 were adopted to validate the model test. The value of the first 10 years is used to calibrate the simulation model, and the second 10 years is used to verify the accuracy of the future prediction. The accuracy of calculation results is relatively high, which denotes that the model is feasible. The R_eff of the observed value and simulated value were 0.84 and 0.80, respectively (Figure 3).

After the model simulation, the simulation results and the efficiency coefficient of each year were obtained. The best R² was 0.84 and the worst year reached 0.80, indicating that the model could explain more than 80% of the flow generation process in the basin (Figure 3). In terms of simulation accuracy, our model was successful in the study area. This benefits from the setting that the model takes into account both the glacier melting runoff in summer and snow water equivalent in winter, which results in the simulated values being close to the actual observed value. However, the degree-day factor method of the glacier melting simulation in the model had insurmountable defects. For example, with the increase of the time resolution, the simulation precision lowered (Hock 2003). Therefore, when comparing the model simulation with the observed monthly runoff depth (Table 4 and Figure 4), the effect was higher than the precision of the daily runoff depth. In addition, the correlation R² between the simulated monthly runoff depth and the observed value was as high as 0.98. In contrast with the annual runoff depth, the error in the maximum annual runoff depth was less than 6% during the model validation (Table 4). Since there is no practical significance for analyzing other months because of the dry season with less runoff, the precision was considered to be acceptable. In summary, the model performs better in the simulation of monthly runoff depth than daily runoff depth, indicating that the changes in the monthly runoff depth obtained by the model were more acceptable.

Besides the influence of precipitation and temperature, the Poiqu/Bhotekoshi basin may be associated with the terrain and water vapor sources crossing the Himalayas, where the summer precipitation is mainly affected by the Indian monsoon (Zhang et al. 1997). As we found in the results, the Poiqul/Bhotekoshi basin had a similar precipitation pattern to the Tamakosi basin, which ensured the accuracy of interpolation using the adjacent meteorological data. The daily precipitation, daily temperature, and daily runoff depth data in the last two decades were simulated and are shown in Figures 5 and 6, respectively.

Rainfall and temperature are the two most important driving factors directly involved in daily operation in the HBV model, so we need to analyze both of them first.

In the past two decades, the daily rainfall in the Poiqu/Bhotekoshi basin changed roughly following the annual cycle. Only in a few months and a particular year was the daily precipitation found to be higher or lower. Some examples are as follows: the daily precipitation in July 2008 reached 140 mm; the whole-year precipitation was higher in 1994 and 1995 (Figure 5); and June–August of every year was basically the abundant precipitation cycle. However, the trend observed from the station data was exhibited differently due to the variations in their altitude. In fact, the precipitation did not always decrease while the altitude got higher. The elevation range from 1,200 to 2,000 m had the precipitation gradually increase as the altitude increased. When the altitude reached a certain threshold, such as 4,000 m or greater, the precipitation was actually very low, and the multi-year average precipitation was only approximately 600 mm. The areas with relatively high precipitation were mainly concentrated on the southern Himalayas, which are primarily covered by coniferous and broad-leaved forest vegetation, and have less arid lands and paddy fields. During July–August, the precipitation was abundant, leading to frequent floods during the rainy season.

Based on the in-situ meteorological observations, kriging interpolation was implemented to map the annual high, annual average, and annual low temperature. After that, a comparison was made among these temperature estimations. The daily temperature changes also followed the annual cycle. The difference between the annual maximum and minimum temperature was 0.22–39.56 °C. However, the temperature increased with the altitude decline. At altitudes between 3,800 and 4,300 m, the annual average temperature was around 3.25–3.82 °C (Figure 6). However, a warming trend was shown on the northern Himalayas, where the temperature tendency rate was 0.32 °C × 10 a^–1. This may contribute to the accelerated snow melting in that region. When the altitude was approximately 2,000 m, the annual average temperature reached 14.56 °C. At altitudes between 1,300 and 1,400 m, the annual average temperature was 17–18.56 °C, which was higher than that at 2,000 m by 4 °C, and its maximum reached 36 °C. Furthermore, the temperature tendency rate was only 0.006 °C × 10 a^–1, and the elevation showed a slowly increasing trend. At the altitude of 900–1,000 m, the average annual temperature was 21.34 °C. At altitudes below 100 m near the equator, the temperature was very high, with the maximum temperature being as high as 39.56 °C, and the minimum temperature 2.35 °C.

The multiple regression method was adopted to obtain the daily average temperature of the Poiqu/Bhotekoshi basin based on the observation data of weather stations. At the same time, the daily precipitation data were also obtained. Through mathematical statistics, it was determined that in the last two decades the basin temperature and the precipitation showed upward trends with average changing rates of 0.15 °C × a^–1 and 0.97 mm × a^–1, respectively.

Snow water equivalent in Poiqu/Bhotekoshi basin from 1988 to 2008 was calculated based on Equation (1) (Figure 7). The snow accumulation showed a significant decrease in the region with snow accumulating for less than 120 days during the 20-year period; the regions with increasing snow accumulation had snowy days of around 120–350 (Sun et al. 2014) and the area of stable snow cover was gradually expanding. According to the simulated snow water equivalent, the expanding regions of snow water equivalent were positively correlated with the snow accumulation days to a great extent. The snow water equivalent was increasing year by year while snow accumulation days were decreasing in Poiqu/Bhotekoshi basin (Figure 7), which may indicate a rising temperature causing a transition from perennial snow to seasonal snow while the snowfall was increased. The snow area was shrinking, while the maximal snow areas first showed a fluctuating increasing trend and a decreasing trend after that. The average accumulated snow showed a fluctuating decreasing trend with an annual rate of –1.00 × 10³ m³ × a^–1. There was a gradual increase in snow water equivalent and the rise started at each fall. The drops of average annual snow water equivalent were 10.48 and 11.73 mm in 1992 and 2005 respectively, while the peaks were 17.59 and 22.38 mm in 1994 and 2006 respectively. From 2006 to 2008, snow water equivalent kept soaring, which was consistent with the previous finding of temperature rise and shrunken snow area. In summary, snow accumulation in most of the basin areas is decreasing year by year; on the other hand, the stable snow area is expanding and perennial snow cover is shrinking, which may be influenced by global warming.

Because the soil water content between the upper and lower layers is highly correlated, the soil water content in the upper layer can be estimated by using soil moisture storage in the lower layers through the HBV model. Soil moisture storage showed a nonlinear increasing trend in its depth (Figure 8), being consistent with the complexity of soil moisture of Poiqu/Bhotekoshi basin in a vertical direction. During the 20-year period, the soil moisture storage gradually increased with a daily rate of between 64.12 and 90.62 mm. The average values of winter–spring and summer–fall are 67 and 85 mm respectively, which demonstrates an obvious difference in water storage between dry and wet seasons (Figure 8). Many studies found that soil moisture storage was highly correlated with climate factors (Hou & Wulan 2006; Zhao et al. 2014). However, the correlation is also affected by the storage depth and observation season. The correlation coefficients between climate factors (temperature and precipitation) and soil moisture storage were calculated using SPSS software. We found that soil moisture storage was significantly affected by temperature as well as precipitation in the wet season. The results also illustrated that the rise of temperature drove permafrost and glacier to melt, which in turn increased the runoff.

The amount of soil water loss can be described by the soil infiltration depth in the HBV model. The soil infiltration depth showed a gradually increasing trend (Figure 8). Most of the high values, ranging from 80 to 150 mm, were observed in July or August with an infiltration rate of 0.08 mm/min. For particular months it reached more than 150 mm depth with an infiltration rate at 0.10 mm/min. The soil infiltration depth has a positive correlation with the soil moisture storage, in other words, both the snowmelt and rainfall have direct influences on the soil infiltration depth.

The basin runoff needs to be considered in two soil conditions, high permeability and low permeability. In the high permeability of the upper layer, the water level descends following the decreasing rainfall and soil infiltration without interception, then we define the storage threshold of the upper layer as UGL (Equation (5)). When the upper groundwater storage is greater than UGL, the runoff is determined by the coefficient K₀. If it is less than UGL, the value will be determined by K₁. In the low permeability of the lower layer, the groundwater level is no longer declining, or decreasing slowly, even if the rainfall and other factors are still dropping. The runoff in the lower layer is determined by the coefficient K₂.

The storage in upper and lower groundwater and the runoff in upper/lower groundwater were simulated individually based on Equations (5)–(7). The simulation results are shown in Figures 9 and 10. The upper groundwater in Poiqu/Bhotekoshi basin is in an aquifer formed by metamorphic rocks, granites and shallow-water sedimentary rock of the Himalayas (Kazuo 1993). Because the aquifer has good water yield property, the regular or seasonal rainfall or snowfall are all likely to result in considerable storage volumes. The annual peak of the upper groundwater storage appeared in mid-August and varied between 140 and 200 mm (Figure 9). Within a year, the upper groundwater storage increased from March, reached the peak in August, and then fell back in early November. The lower layer is formed by metamorphism and weak metamorphic sedimentary rocks from Precambrian to Paleozoic. Although the peak value of storage in the lower layer was also observed in August, the value only ranges from 120 to 180 mm. The seasons having an increasing water concentration were from July to early September, while the other months were lower than 50 mm. Compared to the upper level, the volume and recharging period in lower-level groundwater are significantly less or shorter than the upper level. As the lower layer has a larger elastic storage water rate, the contribution of periodic inland recharge on the storage is effectively reduced. In general, the contribution of the lower layer to groundwater storage is around 12–24% of the upper layer.

The total runoff in Poiqu/Bhotekoshi basin was calculated as the sum of the upper and the lower runoff (Figure 10). The water level of the upper groundwater goes up as the rainfall infiltration recharge increases. The correlation between the upper runoff, the storage and rainfall are quite significant. The simulation results show that the runoff in lower groundwater is twice as high as the upper level; on the other hand, the storage in lower groundwater is slightly lower than the upper level. When the surface water level rises, the underwater pressure increases, the water infiltration increases, and the peak value of groundwater runoff curve will be higher. In the season of heavy rainfall, the infiltration of rainwater also increases, and the peak value of groundwater runoff curve will be higher. One of the reasons is that the water level of the storage in upper groundwater is always greater than the lower, which leads the groundwater to flow downwards to the soil aquifer. Another reason is that part of the water flows into the sea due to the hydraulic gradient.

In this study, the runoff was represented by the runoff depth. The daily runoff in the simulated annual cycle shows a similar overall trend. It increases in early January, peaks in August, and then declines in October. In August 1988, the maximum runoff reached 714 mm, and the maximum of other years ranged from 300 to 470 mm (Figure 4). According to the spatial pattern of the runoff depth in the Poiqu/Bhotekoshi basin, the runoff depth gradually increased from north to south. The runoff depth on the northern slope of the Himalayas was relatively lower (0.97–340 mm), while the southern slope had an average of 200–2,100 mm. From May–September, the runoff depth increased until it peaked in August, and the overall monthly average was up to 574 mm with a minimum average of 327 mm. The periods with lower runoff depth were mainly in the dry season or winter, for example, January–April and October–December, as shown in Figure 5.

The change of permafrost in the river basin contributes to less than 1% of the water balance (Zhang & Yao 2002). Therefore, it is reasonable to assume that the annual average soil moisture content in the study area did not change. In the calculation, the impact of the precipitation gradient was considered in the precipitation (P) (Equation (10)); E_p is the model simulation result; and E_g is the basin surface evaporation taken as the fixed value of 137.23 mm (Zhang & Yao 2002). When calculated with Equation (10), the basin mass balance from 1988 to 2008 had an average of –135.85 mm × a^–1. The simulated basin mass balance had a correlation coefficient of 0.77 (α = 0.01) with the observed mass balance (Figure 11) (Pu et al. 1998; Fujita & Ageta 2000). We concluded that the model could effectively simulate each part of the water balance. The simulation results were reasonable, and the estimated model parameters could reflect a relatively true GMB in the Poiqu/Bhotekoshi basin.

Bhanu et al. (2016) conducted a simple hydrological calculation using the annual calculation and absorption from snow meteorological and discharge data, and analyzed the GMB observation values of Indian Himalaya over the last 40 years. Their results are close to the simulated values in our study, and the general trend is consistent with ours. Etienne et al. (2007) used remote sensing technology to estimate GMB. Although the low value range of GMB estimated by the author is between –850 and –700 mm × a^–1, which is lower than the estimated value in our study, the glacier loss rate between 1999 and 2008 is twice that of 1988–1998, which is consistent with the analysis trend in our study. Ann et al. (2015) analyzed the response of GMB to climate change. Their model ignored the physical causality between input and output, which was not conducive to the analysis of the intermediate hydrological process, and our model overcomes this shortcoming.

The calibrated HBV was implemented in the transnational Poiqu/Bhotekoshi basin. The model first simulated the forming and flowing process of the runoff in the river basin, where the LULC could be classified as alpine tundra zone and a mountain vegetation zone. The glacier area accounts for 29% of the total area of the basin, and the effect of glacier ablation on the basin is very important. Therefore, the accuracy of the model for the glacier ablation simulation greatly affects the accuracy of the model for runoff simulation. The main advantage of the HBV model method is that it needs less data. The value of 1988–1998 was used to calibrate the model, and the value of 1999–2008 was used to verify the accuracy of future prediction. R² was 0.84, indicating that the model is feasible. Then, the model performs using some initial inputs, which were the in-situ monthly observations of temperature and precipitation, to simulate the traffic trails as well as to forecast the possible changes of the runoff in response to the climate change. This model can now be used to perform some monthly mountainous traffic trail forecasting based on the Poiqu/Bhotekoshi basin, which is located in a particular zone transitioning from a higher altitude to a lower altitude.

From the analysis of glacier/snow modules, soil modules and groundwater module, the days of snow accumulation in the Himalayas was decreasing, in particular, the permanent snow area is shrinking with an annual rate of –1.00 × 10³ m³ a^–1 because of global warming. The stable snow area was increasing, which results in a slow rise of snow water equivalent. The annual soil moisture storage was mainly affected by air temperature, and the correlations between soil moisture and temperature in spring, and soil moisture and rainfall in summer were the most significant. The model also simulated the storage and runoff of the most important groundwater response area. The simulation was performed on two soil layers, and then the water mass balance values were calculated. The R_eff of the daily runoff depth simulation after the runoff calibration process was above 0.80. Therefore, the HBV model had a good simulation efficiency on the annual daily runoff depth in the Himalayan cross-river glacial mountains. The water and heat conditions of glaciers in the Himalayas are constantly fluctuating, and the annual material balance is also different. Sometimes the accumulated amount is greater than the melting amount, and a positive balance appears, which is conducive to the development of Himalayan glaciers; otherwise, it will produce negative balance, leading to glacier shrinking.

The results show that the simulated values of the HBV model are reliable and can be used to estimate the GMB of cross-border mountain basins with huge elevation difference, and provide scientific data support for water resources management. However, in this study, hydrological observation data and evapotranspiration data were inadequate, and more abundant model-driven data can further improve the simulation accuracy. In addition, hydrological and underlying surface data can be used to regionalize the parameters and provide scientific data support for hydrological forecasting in the future study.

This study was supported by the National Natural Science Funds of China (grant no. 41971226; 41871357), the ‘One-Three-Five’ Project of Chinese Academy of Sciences (grant no. SDS-135-1708) and the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA19030303).

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