Introduction

Pakistan, once a water surplus country is now categorized as one of the most water stressed countries across the globe. The rapidly increasing population has already reached 225 million and is adding to the stress on water resources with annual per capita water availability declining from 5000 m3 in 1951 to 1100 m3 in 2005 and would likely fall below 800 m3 by the end of 2025 with an anticipated population of 250 million people. The country’s annual water demand is increasing at nearly 10% per year which will cause a yearly volume deficit of around 80 km31. This volume deficit is more than half of the Indus River’s total annual average flow2.

During the last decade the world witnessed a prominent increase in the occurrence and magnitude of extreme events. This increase is expected to get more pronounced in future if the anthropogenic activities are not managed accordingly3. Pakistan had been no exception and has been facing intense floods and heat waves in recent years4,5. In the last 5 decades, the annual average temperature in Pakistan has increased by almost 0.5°C. Occurrence of heatwaves has also increased five folds with the recent ones in year 2015, 2018, and 2019 wherein the temperature crossed 50°C. The highest ever temperature for the month of March since 1901 was also observed in year 2022 with 49.5°C6,7. Temperature projections for future show that there will be a rise up to 4.5°C by the end of twenty-first century8. Annual precipitation in Pakistan is highly varied across regions but shows an increasing trend9. Lutz et al.10 also documents an increase in annual, summer, and winter precipitation trend over the Upper Indus Basin. Overall, the future climate would mostly be dominated by increased temperatures, accelerated melt rates, and contrasting but majorly increasing precipitation in the UIB11. The impact of all these climatic changes are more pronounced for Pakistan because of its minimal social and financial adaptive capacity with high population and dependence on natural environment12.

Glacier loss and changes in streamflow have been observed across all mountainous areas around the world; e.g. the Andes, the Hindukush-Karakoram-Himalaya (HKH), and the Swiss Alps13. However, more research is needed to understand the glacier retreat and dynamics, and reduce uncertainties in mass balance assessments14. The need for understanding glacier retreat and dynamical streamflow change is important for millions of people living in downstream regions, since they rely on meltwater to sustain their livelihoods15,16,17,18,19.

The Himalayas and the Tibetan Plateau host the origin of Asia’s most important rivers and the livelihood of a large part of the continent depend on them. Upper Indus Basin (UIB) is a transboundary river basin, covering parts of China, Pakistan, India, and Afghanistan. The UIB also sustains the world’s largest contiguous irrigation system, about 70–80% of which is provided by snow and glacier melt10. With varying precipitation patterns, the monthly streamflow variation is very high in the basin. The uncertainties associated with the future response of the high elevation HKH region and the large dependency of the water availability in the Indus Basin on these upstream sources of water, makes this basin a climate change hotspot20,21.

In the past decade and presently, the whole Indus basin especially the UIB has been the focal point for research for the climate impact studies and to understand the complex hydro-climatology of the region22,23,24,25,26. There are still; however, uncertainties associated with the results because of varied glacier response in the region. Himalayas exhibit the fastest melt rate in the past decade, while the Karakoram and Pamir ranges show little to no change in glacier mass27.

As per the global climate projections, the temperature in UIB is expected to rise and accelerate the glacier melt28. However, there are contrasting evidence of climate change from different regions of UIB29. The UIB historic temperature records show a strong increase in winter temperature, but a decreasing trend for summer months in some regions30,31. Hussain et al.32 reports increasing trends for seasonal Tmax in winter, spring and autumn at rates of 0.38, 0.35 and 0.05 °C/decade, respectively, while decreasing in summer with −0.14 °C/decade. Also, Tmin exhibits significant decreasing trends in summer and autumn.

Reggiani and Rientjes33 have discussed the hydrological water balance of the Hindukush-Karakoram-Himalaya (HKKH) system shared by Pakistan, India, and China. The streamflow at Besham Qila, just upstream of Tarbela reservoir, has a slightly decreasing trend in contrast to the increased ice melt in the region. It is important to mention that Indus gets 45% of glacial melt from UIB, the highest share among ten rivers of Himalaya. Annual outflows from strong-glaciated sub basins like Hunza, Shyok, and Shigar are decreasing or either stable suggesting decreased precipitation or its storage as snow. However, from the overall trend analysis of 50 year basin outflow records, increased outflow volume is not recognizable.

Most of the impact studies carried out previously were based on conceptual models focusing majorly on simulation of streamflow and extreme events34,35. Lumped or conceptual hydrological models, such as Stanford watershed model36, Snowmelt Runoff Model37, University of British Columbia (UBC) Watershed Model38, Hydrologiska Byrans Vattenbalansavdeling (HBV)39, and others, are widely used for streamflow prediction because of low computational cost, their simplicity and lower data requirements40,41. These models, however, represent the hydrological processes with much less details42, and use averaged parameters which provide rational results only when there are no substantial changes in watershed conditions43.

The glaciers in the Hindukush, Karakoram, and Himalayan mountains are one of the largest of the world. These glaciers provide water constituting more than 80% of the annual river flows, and caters to the needs of people of eight countries from Afghanistan in the west to Myanmar in the east. Pakistan’s reliance on UIB for its food and economic security together with the loss in storage capacities of its large multipurpose storage structures has increased the water challenges in the country. It is; therefore, imperative to study the hydrological regime in the region under the impact of climate change, using improved hydrological modeling tools.

The unique aspect of this study is to contribute to improving the understanding of the hydrological processes in the Upper Indus Basin and assessing their response to changes in climate in terms of any changes in future water resources, using coupled MIKE SHE/MIKE HYDRO RIVER model. MIKE SHE is a fully distributed, physically-based hydrologic modelling system from Danish Hydraulic Institute (DHI) developed by Refsgaard & Storm44. For reliable representation of surface water-river flow interaction and dynamics, the hydrological model was coupled with a hydrodynamic river model, MIKE HYDRO RIVER (a newly proposed version of MIKE 11). MIKE SHE with its distributed and physically-based routines has been widely used globally45,46 and for smaller individual catchments within the UIB47,48. The coupled MIKE SHE/MIKE HYDRO RIVER model would be utilized for first time for the whole UIB in this study to quantify the projected changes in mean annual stream flow and seasonal changes under the influence of climate change, up to the end of twenty-first century.

Methodology

Study area

The UIB is located in the Hindu Kush, Karakoram, Himalaya mountain ranges and on the Tibetan Plateau as shown in (Fig. 1). The figure was generated using QGIS v3.34.4 (https://www.qgis.org/download/). Seven tributaries of the Indus River drain through the UIB, covering a drainage area of ~ 170,000 km2 at the Tarbela reservoir in Pakistan. The UIB has a maximum altitudinal extent of 8500 m, with average elevation of 3750 m. UIB is one of the most glaciated places on the face of Earth, with 22,000 km2 of glaciated area49,50.

Fig. 1
figure 1

The Upper Indus Basin draining at Tarbela reservoir in Pakistan.

Full size image

The UIB is an important water resource of about 170,220 million m3 of water annually51. The annual snow and glacier melt have significant contribution of up to 60–70% to the water availability every year52. The precipitation patterns in the basin effects the hydrology and overall water resources for most parts of Pakistan. Climate induced changes in precipitation patterns can,therefore, severely affect the water availability in Pakistan53,54,55.

The climate pattern of UIB is complex especially due to diverse topography in its different parts. The weather and climate are influenced by the summer monsoon system and western disturbances. The valley regions of the basin are mostly arid, with mean annual precipitation of 100–200 mm. The precipitation; however, increases with altitude and can rise up to 600 mm at 4400 m elevation, with glacier accumulation amounting to around 2000 mm above 5000 m altitude56. Much of the annual precipitation is during spring and winter months due to western circulations57. The local circulations and monsoon system accounts for one third of the total annual precipitation. The monthly temperatures show that July is the hottest and months of December-January–February are the coldest58.

UIB discharges an average of 2200 m3/s of water annually just downstream of Diamer Bhasha dam at Shatial Br. gauge station. Despite the various claims of increased glacier melt in the region, the flow data of past 30 years show consistent and near stable trend in the basin, as shown in (Fig. 2). Most of the water available in the basin is in the monsoon season, from June- August. The month of July receives the major share of annual water in the basin with an average flow of approximately 6500 m3/s.

Fig. 2
figure 2

Mean annual flow at Shatial Br. Station.

Full size image

Data collection

The data required for the study were divided into physiographic and climate categories. The physiographic data includes digital elevation model (DEM), Land cover data, and soil datasets. These datasets have been acquired from various web sources of relevant departments and formatted into the required model format. More detailed information on these datasets is covered in section 2.4.

For UIB, the meteorological data for five stations were used which include daily precipitation and maximum and minimum air temperature. The details of stations are presented in (Table 1).

Table 1 Meteorological Stations in the study area.
Full size table

The observed discharge data used for the calibration and validation of hydrological model was acquired from Water and Power Development Authority (WAPDA). The dataset was for Shatial Bridge gauge station, recorded at daily time step, installed just downstream of Diamer-Bhasha Dam for the period 1984–2013.

To assess future water availability in the region, the hydrological model requires future climatic datasets. The future climatic data has been obtained from CORDEX-SA (0.44°) simulations. These datasets were obtained under the CORDEX-SA project by dynamically downscaling CMIP5 GCMs using various RCMs. For this paper, these datasets were accessed from earth system grid federation (ESGF) node at https://esgf-data.dkrz.de/projects/esgf-dkrz/ (ESGF-DKRZ 2014). Table 2 lists down future climatic datasets used in this study. The GCMs were selected to cover the broad range of climate predictions for HKH region which include warm and dry climate, warm and wet, cold and dry, and cold and wet climate for future years. The downscaling RCM was chosen to be same to consider changes across GCMs only.

Table 2 Driving GCMs and downscaling RCM model.
Full size table

The well-tested CORDEX-SA-downscaled data of CMIP5 GCMs was used because of its wide applicability in the UIB. Moreover, the study mainly aims to assess response of the UIB hydrological system to any changes in climate using MIKE SHE, and not the efficiency of CMIP6 or CMIP5 GCMs in simulating UIB’s climate. Some other recent studies have also employed CMIP5 datasets like59,60

Model design

With physically-based distributed models, one can work around with many of the deficiencies inherent to use of lumped-conceptual models. This is primarily done by usage of parameters having a physical relationship and can include spatial and temporal variability to the parameters61 to a greater extent as compared to lumped models. Therefore, this study has also selected MIKESHE, which is a physically based and a fully distributed model.

MIKESHE is a powerful, physically based, fully-distributed model for 3D simulation of watershed scale hydrologic systems. It models entire land-based flow of water within the watershed. The model has been used globally for studies ranging from small watersheds to large river basins for regional studies62,63,64,65,66. Despite of the vast applicability of MIKE SHE to varying topographical, and hydrological regimes, it has never been utilized for hydrological simulation of the highly complex and important river basin, the UIB.

Sahoo et al.43 discussed systematic procedure for the calibration and validation of MIKE SHE for Manoa–Palolo mountainous drainage basin on the island of Oahu, Hawaii. The model produced reliable simulations of streamflow and groundwater levels. Sultana and Coulibaly67 evaluated the suitability of a physical based, fully distributed coupled model MIKE SHE/MIKE 11 to predict changes in hydrology of Spencer Creek watershed, Ontario, Canada. The results suggest model adequacy for climate impact studies at watershed scale.

Jaber and Shukla68 believe that the coupled MIKE SHE/MIKE 11 model can simulate complex watershed processes much more efficiently than many other models. The MIKE SHE/MIKE 11 (recently improved to MIKE HYDRO RIVER) coupled model had also been selected by various other researchers because of its capability to model the dynamic exchanges between the surface flow, aquifer system, and the river system. Refs.69,70,71.

MIKE-SHE is based on finite difference approach for simulating overland and river flow, and saturated (groundwater) and unsaturated (sub-surface) zone flows, employing diffusion-wave approximation. The study area is discretized into a network of grids. The resolution of these grids depends on the user input based mostly on the available computation power. For this study, the area was discretized into 1000m X 1000m grid due to limited computational power available for these simulations. The evapotranspiration, interception storages and snowmelt routines are based on analytical solutions72. The coupled MIKE-SHE/MIKE 11 model allows using the complete Saint–Venant equations73. During simulation, water levels within coupled streams are transferred from MIKE-11 to MIKE-SHE river links, at each computational step. MIKE-SHE estimates the river-aquifer exchange and overland flow for each river link and fed them back to MIKE 11 as lateral inflow or outflow72. MIKE-She’s robust interface also allows user to select appropriate computational time step and input and output data time steps74.

Similar to many other watershed models, the MIKE SHE model uses degree-day approach and precipitation and temperature lapse rates (to account for change in precipitation and temperature along the elevation of the watershed) to model distribution of these weather parameters across the watershed and to calculate snowmelt. In this study, the snowmelt module has been used to calculate glacial melt by using the assumption (similar to other published studies e.g. Hashmi et. al.75) that the melt water received during July–September is mostly glacier-melt., However, the model lacks any specific glacier module.

Modeling processes

Model application procedures for MIKE SHE includes time series and raster data development, setting up saturated and unsaturated zones, coupling with MIKE HYDRO RIVER for the simulations of river dynamics, and model calibration and validation.

The daily precipitation and temperature time series data for the years 1999–2008 were edited within the time-series editor module of MIKE SHE. The station-based point rainfall and temperature data were spatially distributed over the entire basin using area weighted time series method in MIKE SHE. The daily time series of Potential Evapotranspiration (PET) was calculated by a temperature‐based equation76 using the daily average temperature (derived using daily minimum and maximum temperature time series) outside the model domain due to limitation of meteorological data.

Digital elevation model (DEM) is the 3D representation of terrain or topography of the area. The DEM for the watershed was obtained from shuttle radar topographic mission (SRTM) of 1 arc-second resolution (30m). Automatic Delineation processes under ArcHydro extension in ArcGIS was used to delineate DEM grid. The tools create two separate layers i.e. stream network and sub-basins. These GIS layers were manually edited and formatted based on topographical analysis of region from Google Earth, and literature review to account for the natural features overlooked in delineation tool such as closed basins of Tibetan Plateau. The delineated stream network was imported in MIKE HYDRO RIVER interface. Automatic Cross-section creation tool had been used to create cross-sections for the 22 delineated streams for the study area.

Land use data is used to define the hydrologic parameters associated with each land category within the study domain. The land use layer used in model setup is GlobCover (2009). The layer was processed under ArcGIS routine and transformed to MIKE SHE acceptable format. The soil cover was obtained from digital soil map of world (DSMW) developed by the United Nations food and agriculture organization (UN-FAO). Sandy clay loam is the most dominant soil type in the region. The saturated and un-saturated zones have been classified based on the depth to water table, details of which have been obtained as a result of a personal communication with a senior researcher in the field of hydrological modelling analysis working in the UIB area.

Calibration & validation

Successful development of required data files was followed by simulation with MIKE SHE/MIKE HYDRO RIVER coupled hydrologic and fully hydrodynamic river model. The model parameters were evaluated and adjusted based on observed daily streamflow Defining values of various parameters helps in proper demonstration of watershed behavior and its hydrological cycle. The parameters have been tested for their sensitivity to model outputs and calibrated thereafter to achieve a good correlation between observed and simulated results. To evaluate the parameters having most influence on model outputs, sensitivity analysis was done. Sensitivity analysis helps to analyze the realistic parameter ranges and their relationship amongst each other77.

The sensitivity analysis was carried out using both automatic and manual calibration procedures. The identified sensitive parameters and their descriptions for MIKE SHE are tabulated in (Table 3). The adjustment of these most sensitive parameters was followed by calibration and validation of the coupled model. The model has been calibrated and validated for years 1998–2001 and 2002–2004 respectively. To assess the reliability of calibrated model to be used for impact assessment, two different matrices, Nash–Sutcliffe efficiency (NSE) (Eq. 1), and coefficient of determination (R2) (Eq. 2) were used in this study. These matrices with their inherent limitations are still most widely used to predict accuracy of model results. The study has also used kling-gupta efficiency (KGE) (Gupta et al.78). The KGE ranges from − ∞ to 1, a performance above 0.75 and 0.5 is considered to be as good and intermediate, respectively (Eq. 3). The ideal coefficient value is one.

$$NSE = 1 - \frac{{\mathop \sum \nolimits_{i = 1}^{n} \left( {O_{i} - P_{i} } \right)^{2} }}{{\mathop \sum \nolimits_{i = 1}^{n} \left( {O_{i} - \overline{O}} \right)^{2} }}$$
(1)
Table 3 Most sensitive calibration parameters.
Full size table
$${R}^{2}= {\left[\frac{\sum_{i=1}^{n}({O}_{i}- \overline{O })({P}_{i}-\overline{P })}{{({\sum_{i=1}^{n}({O}_{i}-\overline{O })}^{2})}^{0.5}{({\sum_{i=1}^{n}({P}_{i}-\overline{P })}^{2})}^{0.5}}\right]}^{2}$$
(2)
$$KGE=1-\sqrt{{(r-1)}^{2}+{(\alpha -1)}^{2}+{(\beta -1)}^{2}}$$
(3)

Since the annual flows in the UIB depend more than 50% on glacial and snow melt, the snow melt routine was very sensitive and crucial in model calibration. MIKE SHE employs a degree day model (DDM) for snow melt prediction. To account for seasonal variation in snow melt, a time varying degree day factor (DDF) was used. The DDF starts from 1.5 mm/day per °C in April and is ramped up to 5 mm/day per °C in July, decreasing up till 0.5 mm/day per °C in September. The calibrated values of DDF are similar to those reported by Fatima et al.79, and Latif et al.80. The other sensitive parameters include temperature and precipitation lapse rates to account for the contrasting elevation zones in the region and orographic precipitation. The initial value for Temperature lapse rate was derived from Mukhopadhyay & Khan81, wherein values ranging from 0.5 to 1°C/100m was reported for different sub-basins of UIB. For precipitation lapse rate, estimates based on orographic correction factors derived by Khan and Koch82 was referred and manually calibrated for this study. The difference in values may be attributed to lower data representation in this study. Base flow time constant and Manning’s M were also sensitive especially at monthly time scales. These were adjusted to account for flow recession in winter and early spring months, and surface characteristics, respectively.

Future climate data generation using empirical quantile mapping

Future climatic data time series to be input to the validated hydrological model was generated by bias-correcting the CORDEC-SA climate model data against the observed datasets using the empirical quantile mapping (EQM) method83. This analysis particularly focuses on the mean flow changes, on annual and seasonal scale, up to the end of the twenty-first century under all climatic scenarios. The EQM method is based on empirical cumulative distribution functions built on daily time scale (ecdfs). In this method (Eqs. 4 and Eq. 5), a quantile of the present day simulated distribution is replaced by the same quantile of the present-day observed distribution. The main assumption here was that the bias correction is the same for the same probability function, both in the control period and in the future scenarios. The equation for corrected temperature and precipitation can be expressed as:

$${p}_{hst,m}^{corr}={ecdf}_{obs,m}^{-1}({ecdf}_{hst, m}\left({P}_{hst,m,d}\right))$$
(4)
$${T}_{hst,m}^{corr}={ecdf}_{obs,m}^{-1}({ecdf}_{hst, m}\left({T}_{hst,m,d}\right))$$
(5)

where ecdf−1 is the inverse representation of ecdf.

Future flow scenarios simulation

After successful validation of MIKE SHE model for UIB, the bias-corrected (Sect.  2.6) future climate time series was used as an input to the validated model. The hydrological model was run from 2020 to 2099 to understand the impact of changing climate on regional stream flows.

The future meteorological data was obtained from four different GCMs downscaled by one distinct RCM to evaluate climate predictions across GCMs alone. The GCM datasets were dynamically downscaled to a relatively high resolution ~ 50km to represent the region’s characteristics using Rossby Centre regional atmospheric model, RCA4. The GCMs have been selected after rigorous literature review to have contrasting future climate prediction across the basin. The details of GCMs are presented in (Table 2).

For future climate scenarios, twenty-first century climate projections for an extreme emissions scenario-RCP 8.5 and a medium emissions scenario-RCP 4.5 were used for the period 2010–2099.

With radiative forcing of 8.5 W/m2, the RCP 8.5 scenario is an upper bracket of climate projections. While, RCP 4.5 is a stabilization scenario projecting an increase in temperature for up to 2.4°C by the end of the century and radiative forcing up to 4.5 W/m284

Results & discussion

MIKE SHE/MIKE HYDRO RIVER calibration and validation

Calibration was performed on daily and monthly time steps from 1998 to 2001. The details of calibrated values for identified sensitive parameters are given in (Table 4).

Table 4 Sensitive parameters and their calibrated values.
Full size table

Figure 3 shows calibrated results of the model compared to observed flows. The simulated peaks show a good correlation with the observed data. NS-efficiency and R2 of 0.77 and 0.84, and KGE value of 0.657, respectively, are the evidences that the model results are reliable for further analysis.

Fig. 3
figure 3

Calibration based on daily streamflow comparison.

Full size image

However, the model under predicts the flood peaks in the monsoon season and has slightly over predicted the flows in late autumn i.e. Sep-Oct. This underestimation maybe attributed to underestimated precipitation due to interpolation within the catchment, and the model might not be able to reproduce the glacier melt.

The calibrated model was run for the validation period of 2002 to 2004 without altering model parameters as a standard practice. Figure 4 shows a close match between observed and simulated peaks for the validation period. The model results were acceptable for the validation period with an NS efficiency, R2 and KGE of 0.73, 0.83, and 0.587, respectively. The model captured the seasonal variation in runoff in the Indus River at Shatial Bridge. Khan et al.85 also reports similar calibration and validation results at Shatial gauge station.

Fig. 4
figure 4

Validation results based on daily streamflow comparison.

Full size image

Future water availability

Streamflow from the UIB has a strong correlation with observed temperature. Similar previous studies suggest that this is mostly due to dominance of temperature-sensitive snow and glacial melt in the total flow.

The impact of climate variability and change on future water availability from the UIB was studied using calibrated & validated MIKE SHE/MIKE HYDRO RIVER model under different climate change scenarios, as discussed in Sect. 2.2. Due to limitation of MIKESHE to include a varying snow extent and depth, the modelling exercise assumes that the spatial coverage of snow is constant, even for the future years, when warmer climate is projected across all the climate model projections, and only the amount of snow melt changes with any changes in the temperature.

Inter-annual changes

The ensemble generated in this study projects increased precipitation and warmer temperatures for both RCP 4.5 and RCP 8.5 scenarios. Due to the rise in temperature, the precipitation would fall more as rain than snow, and would also presumably result in increased snow melt contribution.

The discharges simulated by the ensemble run are shown in (Figs. 5, 6). The ensemble flows for both RCP 4.5 & 8.5 show considerable increase (86% for RCP4.5 and 97% for RCP8.5) in mean annual flows when compared with the historic dataset. The RCP 8.5 shows consistently higher predicted flows than the RCP 4.5 because of increased temperatures leading to accelerated snow melt. The ensemble spread also shows more intra-annual variability in flows for the RCP8.5 scenario leading to increased uncertainty. The most prominent difference is in the late century period (2076–2099) where ensemble mean annual flows deviate from historic mean of 2460 m3/s by 101 and 112% for RCP4.5 and RCP8.5 scenarios, respectively. The linear trend lines in Fig. 5 and Fig. 6 also show that the increase in stream flows is projected to be very sharp in the mid-century (2020–2050) compared to the late century where the trend is increasing but flattens. This is probably due to loss of snow cover available in the latter half of century due to accelerated melt in earlier years.

Fig. 5
figure 5

Ensemble flows and spread for RCP 4.5.

Full size image
Fig. 6
figure 6

Ensemble flows and spread for RCP 8.5.

Full size image

The analysis of stream flows for 1st percentile and 99th percentile was done to compute future High daily flows and low daily flows compared to the reference or historic flow at Shatial gauge station for the period 1984–2013 (Table 5). The results show consistently increasing trends for all prediction models across both scenarios. The most prominent increase was for GFDL ESM2M which predicted significant increase in mean, low daily, and high daily flows. The percent change was most evident for low daily flows across all models and scenarios due to increase in warmer temperatures leading to increased snow melt. The maximum increase was reported for low daily flows where GFDL model predicts a 500% increase from reference. The mean daily flows though increasing, show a slightly less percentage in far future for CanESM and IPSL model for RCP 8.5. The results, however, are different than what is reported in other studies (Khan et al., 2020). This raises questions on the uncertainty associated with future models output and also the bias correction techniques applied.

Table 5 Mean daily flows, 1st percentile (low flows), and 99th percentile (high flows) for different model runs.
Full size table

Intra-annual changes

All models, across both the RCP scenarios, analyzed for this study, projected higher long-term mean monthly streamflow compared to the baseline period for the UIB. The monthly discharges simulated under ensemble 4.5 and ensemble 8.5 are shown in (Fig. 7). The analysis of mean monthly river flows shows that the flow starts increasing from mid of March as temperatures increase, in contrast to the reference where flow usually increases from mid-April. The maximum flows are observed near the end of July and decreases to end of year. The average monthly flows have increased by more than 300% in the winter months. A sharper increase for both RCP 4.5 and 8.5 during May–June is due to increased snow melt contribution due to higher temperatures.

Fig. 7
figure 7

Ensemble monthly flows for RCP4.5 and RCP8.5.

Full size image

Assumptions and limitations

Modeling future hydrology for a complex and highly glaciated river basin, such as UIB, is a difficult task, limited by uncertainties associated with both in-sufficient or inefficient quality and coverage of observed datasets, watershed model structure and future climate projections. For calibration and validation of MIKE SHE model, five meteorological stations were used (Table 1) for a short duration due to limited data availability. These datasets are not enough to fully represent the complex terrain (varying elevation zones) and cryosphere processes of a large area such as the UIB. There are no weather gauging stations available in regions of high altitude. Also, a longer span of observed data can further improve the calibration process due to representation of various high and low flow years. For future predictions, the well-tested CORDEX-SA-downscaled data of CMIP5 GCMs were used because of their wide applicability in the UIB, and also because the study focuses on suitability of MIKE SHE for impact studies within UIB.

The coupled MIKE SHE/MIKE HYDRO RIVER adequately models most of the hydrological parameters as well as captured the seasonal variation, during the calibration and validation periods. The stream flow magnitude is; however, underestimated in summer months which can be attributed to relatively inadequate representation of glacier melt processes. The glacier accumulation and melt are not only different from that of snow but these are also different for debris-covered glaciers (which are abundant in HKH region) and clean-ice glaciers. The different types of glaciers in the region and uncertainty about the depths of glaciers86, makes glacier modelling complex and is; therefore, beyond the scope of this study. The misrepresentation of UIB’s cryosphere and climatic processes due to unavailability of in-situ datasets has also compromised the impact study results.

Conclusion & recommendations

Conclusion

The current study has tried to investigate the efficacy of using MIKE SHE/MIKE HYDRO RIVER coupled model in replicating the hydrological processes of the UIB and to assess impact of future climatic changes projected by a multi-model ensemble for two future scenarios up to the end of twenty-first century. The future projected precipitation and temperature in the UIB shows an increasing trend across all scenarios. The rainfall and runoff show a consistent increasing trend for both RCP 4.5 and 8.5 scenarios for all of the UIB, which is relatively more prominent (112% increase from historic mean) for the RCP 8.5, and near the end of twenty-first century. The increase in annual streamflow for multi-model ensemble is 86 and 97% for RCP4.5 and RCP8.5 in response to an annual average warming of 2.3 and 4°C for these scenarios, respectively. This is believed to be due to increase in snow melt contribution due to higher temperatures. It can also be due to declining snow mass in region due to increased precipitation in the form of rainfall and accelerated melt rates. Base flow, overall discharge, and the snow melt seem to increase throughout the next 80 years.

The low and high flows are projected to increase drastically, consistent in all scenarios and models. The historic average would be now comparable to future lows i.e. the otherwise average annual flow would now most likely be the minimum flow in the basin. This indicates an abundant or rather increased water availability and supply up to the end of twenty-first century. This also increases the chances of future floods by many folds. However, the other dimension of this increased water availability raises concerns over the available supplies in the context of loss of much of glacial cover by then. Considering glacier melt, the future water availability in the UIB will increase until substantial decline in the glacial and snow mass.

Pakistan is very susceptible to climate change impacts. It is therefore, imperative to be able to predict climate change impacts and take relevant policy decisions to mitigate the effect of these changes for country’s social and economic well-being. Pakistan’s current annual water demand is 165 Billion Cubic Meters (BCM), with a total water availability equivalent to 290 BCM which is to be exhausted by domestic, agricultural, and industrial demands. Pakistan currently stands with a water surplus however, has an inadequate storage and management system due to which the surplus results in floods or drains down to sea. Increased siltation of Tarbela and Mangla reservoirs has reduced their respective storage capacity by 20–30%. The current storage capacity of Pakistan stands at 15 BCM and stands for 30 days only, and is essentially in-sufficient to help in drought years. This will be reduced further due to loss in live storage of Tarbela and Mangla. This in turn would result in reduction in water availability for agriculture especially for growing seasons of Rabi and early Kharif, and reduced hydro power energy potential87.

With projected annual average river flow of 4500 m3/s (equivalent to volume of 142 BCM from UIB alone) under RCP8.5, and historical mean annual flow of 2200 m3/s (69 BCM) at Shatial gauge station, the study concludes that there’s enough water available in the Indus Basin to cater for the needs of Pakistan for the next 80 years, provided the country invests in water infrastructure. In this context, the under construction Bhasha and Mohmand Dam seems promising steps for Pakistan’s growth, adding up to 12 BCM of additional storage by the end of 2030. For an agro-economic country like Pakistan, timely water distribution is important for the crop-productivity in the vast Punjab and Sindh plains. Bhasha dam’s geographical location and capacity will ensure more control over the volume and timing of distribution of water. The Bhasha dam is also believed to ensure extended service life for Tarbela dam and Dasu Power project by capturing much of the incoming silt load.

Recommendations

The increased availability of water predicted during the summer monsoon can cause catastrophic floods or loss of valuable water to the sea if not stored. This study provides some valuable insights on importance of storage structures, both upstream and downstream, to cater the uncertainty of water resources in the UIB. Large dams are believed to be the best adaptation strategy that could address water, energy, and food security simultaneously.

Despite the above findings, the insufficient understanding of the glacier water balance, together with the complex climatic conditions and the varying topography, make future predictions on water resources availability difficult in the UIB. The breadth and complexity of hydrological modelling, particularly for UIB, cautions stake holders that sometimes even a well calibrated model could end with ambiguous inferences.

Therefore, more research is needed to evaluate the uncertainties in approach used. The uncertainty in model predictions and observed datasets needs further evaluation, with other RCMs and hydrological models preferably. Future research should also be towards implementing an ice/glacier routine into physically based hydrological models, which is a limitation of MIKE SHE used in present study. Also, the spatial coverage of hydro-meteorological data grid is limited as the gauging stations are very few compared to the requirement. This is also true for surface water, climatic data, and groundwater levels. The inadequate and ill planned gauging system hinders reliable assessment of the available resources. In-situ monitoring stations should be established particularly for glacier monitoring.