Modelling climate change impact on water resources of the Upper Indus Basin

The Upper Indus Basin (UIB), situated in the Hindukush Karakoram Himalaya (HKH) region, supplies water to one of the largest irrigation systems in the world, meaning that climate change threatens the water, food and hydro-energy security of 200 million people through interruption of water resources (Immerzeel et al. 2010; Amin et al. 2018; Scott et al. 2019). Furthermore, people living within the UIB are vulnerable to floods during the monsoon period, which cause significant economic losses and casualties each year (Shrestha et al. 2015). Any climate change-driven alteration to the intensity, spatial extent or timing of the monsoon, therefore, will also have significant impacts on people in the region. Hence, it is crucial to be able to predict future changes in climate patterns and their impact on river flows of the Indus River at a regional level for sustainable water resource management (Ali et al. 2009; Clifton et al. 2018).

The global average temperature is estimated to increase by 2.1–3.5 °C by 2100 in the intermediate greenhouse gas emissions scenario (SSP2-4.5) (IPCC 2021), which will have a significant impact on global precipitation resulting in more extreme patterns of longer heavy rainfall events followed by long dry intervals (Zhang et al. 2013). In future, warming in the HKH region is predicted to be at least 0.3 °C to 0.7 °C higher than global warming if global warming is kept at 1.5 °C (Krishnan et al. 2019). The amount of precipitation and the timing of river flows in the region are a function of two climatic systems, the Indian monsoon system and the winter westerlies disturbances (WWD). The hydrometeorology of the western Himalaya, including the UIB, is dominated by heavy winter snowfalls from the WWD, and thus, the major contribution to river flow comes from snow and ice melt. Therefore, most studies pertaining to the present and future water resources of the Indus River focus particularly on snow and glacier change in the Karakoram (e.g. Immerzeel et al. 2012; Bocchiola & Diolaiuti 2013; Lutz & Immerzeel 2013; Ragettli et al. 2013; Soncini et al. 2015; Garee et al. 2017).

The results of previous studies from the UIB show inconsistency with respect to both the hydrometeorology (e.g. Fowler & Archer 2006; Khattak et al. 2011; Hasson et al. 2016; Bashir et al. 2017) and predicted future changes (Hasson et al. 2014; Lutz et al. 2016), which complicates sustainable water management efforts. For instance, an increase (Mahmood & Jia 2016) and a decrease (Babur et al. 2016; Lutz et al. 2016) in water yield were predicted during summers in the Upper Jhelum River Basin, while Fowler & Archer (2006) found striking differences in temperature behaviour of the UIB. A consistent increase in the diurnal temperature range was also observed in the UIB (Forsythe et al. 2014). Khan et al. (2019) reported a faster increase in the minimum temperature (0.17–0.37 °C/decade) than the maximum temperature, while Iqbal et al. (2016) observed a faster increase in the maximum temperature than the minimum temperature in the UIB. Differences in temperature and precipitation projections lead to significant variation in discharge estimates from the UIB (Akhtar et al. 2009; Lutz et al. 2016). Predicting the future water supply from the highly glaciered UIB is a challenging task due to highly uncertain climate projections, low-quality reference climatic data, complex climate systems and topography, and model uncertainty (Khan et al. 2020).

In this study, the Soil & Water Assessment Tool (SWAT) is used with newly reanalysed datasets at finer spatial and altitudinal representations than previously available, and with the selection of the most suitable future climate models for this region. Our approach is to look specifically at variations in projected temperature, precipitation and river flow at the sub-basin scale, rather than at the scale of the entire UIB basin, something that has been less reported on in previous studies. We examine spatial and temporal changes in water balance components including water yield, evapotranspiration and soil moisture in the UIB; this work enhances understanding of hydrological processes over the UIB. We examine the predicted changes in temperature and precipitation from the global climate model (GCM) and then model the impact of those changes on water resources by using GCM data to drive a hydrological model. More specifically, this study assesses the effects of climate change on water balance components by parameterizing the SWAT model and evaluating its performance in specific glaciered and non-glaciered sub-basins of the UIB through comparison of outputs against measured flow data. The warming implications for snow/glacier melting from the sub-basins of the UIB are assessed by comparing against baseline and historical data. Simulations are made for the near future (2020–2050) and by the end of the century (2070–2100) under two different RCPs.

The Indus River originates on the Tibetan Plateau, flows through India and enters Pakistan where it provides most of the country's irrigation water (Charles 2015). This study focussed on the UIB which includes the western Himalaya, Hindu Kush and Karakoram mountain ranges. Within the UIB, three sub-basins were identified where the river is gauged and which vary in terms of their climate and topographic characteristics (Table 1): Upper Indus River Basin (UIRB) at Tarbela, Upper Jhelum River Basin (UJRB) at Mangla and Kabul River Basin (KRB) at Nowshera (Figure 1). The Kabul River originates in the Hindu Kush Mountains and receives much of its flow from monsoon precipitation. The Jhelum River flows from the Himalayan ranges and monsoon precipitation, which is at its highest during the summer, contributes significantly to river flow (Rizwan et al. 2019). The Indus River itself receives most precipitation during the winter from the western disturbances but little in the summer as almost 90% of its catchment is in the rain shadow of the Himalayan belt (Immerzeel et al. 2009). In addition, it is fed by major tributaries, such as the Shyok and Hunza rivers, whose flow derives mainly from the melting of snow and glaciers (Immerzeel et al. 2015).

The monsoon precipitation influence is greatest along the Himalayan arc and decreases towards a north-western direction where mid-latitude westerlies have a strong impact on the climate of north-western region of the UIB (Palazzi et al. 2013). WWDs contribute mainly to winter snow (Hewitt 2011; Ali et al. 2015) that melts in spring and contribute to the flows of the Indus and Kabul rivers.

River discharge data were obtained from the Water and Power Development Authority (WAPDA) and Pakistan Meteorological Department (PMD). For stations in the UIB at Tarbela, Mangla and Nowshera, monthly discharge data for the period of 1990–2013 were acquired. The mean monthly runoff of UIRB, UJRB and KRB varied between 1,039–6,065, 263–1,806 and 205–1,596 m³ s⁻¹, respectively, for the study period. In the UIRB and KRB, the highest flows usually occur in June and July, whereas the lowest flows occur in December and January; in UJRB, the highest flows occurred a little earlier in May, likely driven by the spring melting of winter snow.

The SRTM DEM (Shuttle Radar Topography Mission Digital Elevation Model), with a spatial resolution of 90 m, was used for delineation of watersheds of the UIB using ArcSWAT in the ArcMAP10.5 environment. ArcSWAT is freely available, and ArcSWAT extensions can be downloaded from the US Department of Agriculture–Agricultural Research Service website (https://swat.tamu.edu/software/arcswat/). Land-use/land-cover (LULC) data were obtained from the European Space Agency (ESA). The ESA-GlobCover data comprises 23 LULC classes, available annually from 1990 to 2015 (Arino et al. 2007) (Supplementary Material, Figure S1c).

Soil data were obtained from the Harmonized World Soil Database (available from http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/) (Supplementary Material, Figure S1a).

GCMs are the primary source of information in assessments of climate change impacts at regional to global scales (Gebrechorkos et al. 2019). The subject of greater concern is that climate models are still unable to reproduce precisely relevant climate processes (Knutti & Sedláček 2012). However, all climate models are consistent and considered to be more accurate in projecting temperatures, whereas precipitation projections are highly variable (Lutz et al. 2016). Outputs from GCMs can show a high degree of variation, leading to uncertain predictions of future hydro-meteorological conditions (Lutz et al. 2016). Most GCMs poorly represent the South-Asian monsoon regime (Turner & Annamalai 2012). The raw climate outputs from GCMs are rarely used because of biases due to limited spatial resolution, incomplete knowledge of climate system processes and simplified physics and thermodynamic processes (Mehrotra & Sharma 2019). Biases or errors in raw climate model outputs can also be due to historical observations (Ramirez-villegas et al. 2013). Hence, it is important to bias correct GCM simulations in order to produce reliable future climate projections.

Future monthly streamflow was simulated with five GCMs: GFDL-ESM2M designated as GCM1, HadGEM2-ES as GCM2, IPSL-CM5A as GCM3, MIROC as GCM4 and NoerESM1-M as GCM5, along with two emission scenarios (RCP4.5 and RCP8.5) for future projections (2030–2050 and 2070–2099) relative to the historical period (1974–2004) (Table 2). Sperber et al. (2013) explained that no single model can satisfactorily represent the monsoon regime over the HKH region. Results from individual climate models (GCM1, GCM2, GCM3, GCM4 and GCM5) were combined together as a multi-model mean (MMM) to provide a more robust estimate of the change in climatic variables and water balance components. A multi-model ensemble of statistically downscaled GCMs over the Indus River Basin was used, as it was likely to capture precipitation and temperature climatology more effectively than individual GCMs (Su et al. 2016).

In this study, the Climate Forecast System Reanalysis (CFSR) data of minimum temperature, maximum temperature and precipitation over the UIB from January 1980 to December 2013 at a monthly scale were used for constructing the GCM downscaling models to produce bias-corrected temperature and precipitation. The downscaled GCM data of temperature and precipitation were bias-corrected against CFSR data acquired between 1980 and 2013. CFSR data were chosen for downscaling and bias correction of GCMs due to their relatively fine spatial resolution (∼38 km²) in line with study requirements. The climate change toolkit (CCT) was used to extract, downscale and bias correct the outputs from GCMs (Vaghefi et al. 2017). The Bias Correction Statistical Downscaling (BCSD) module is provided in the CCT which is used to bias correct GCM data against CFSR data (Supplementary Material, Figure S9).

Hydrological models are used to study climate change impacts on water resources at global (Van Griensven et al. 2012; Bressiani et al. 2015), regional (Garee et al. 2017; Adnan et al. 2019), and watershed and ecosystem scales (Bhatu & Rank 2017; Terskii et al. 2019). Many models have been developed that simulate river flow and runoff, from small catchments, right up to global scale (Sitterrson et al. 2017). Models to help assess and predict hydrological changes include the Hydrologic Simulation Program Fortran (HSFP; US EPA 2019), the SWAT (Arnold et al. 1998) and the Hydrologic Engineering Centre Hydrologic Modelling System (HEC-HMS 2000). Here, the SWAT model was used, as it includes a component suitable for modelling snowmelt hydrology in watersheds where river runoff is strongly correlated with spring snowmelt (Pradhanang et al. 2011; Abbas et al. 2019), such as in Alpine catchments (Grusson et al. 2015; Tuo et al. 2018), and has efficiently simulated streamflow in the UIB previously (e.g. Babur et al. 2016; Garee et al. 2017; Shah et al. 2020).

The SWAT model functions on a continuous time step and was developed by the U.S. Department of Agriculture Agriculture Research Service (USDA-ARS) (Arnold et al. 2012). The model uses the soil conservation service runoff curve number (CN) method and degree-day factor to compute surface and snowmelt runoff, respectively (Duan et al. 2018a, 2018b). Model components include agriculture management, hydrology, weather, plant growth, erosion/sedimentation, pesticides, nutrients, channel routing and pond/reservoir routing (Arnold et al. 1998). It is a physically based, semi-distributed model that requires primary input data on topography, weather, soil properties, and vegetation and land management practices in the watershed. Daily potential evapotranspiration is computed using a temperature index method (Hargreaves & Samani 1985).

The SWAT model has been used for hydrological modelling studies across the globe (Duan et al. 2018a, 2018b; Tuo et al. 2018; Tan et al. 2019). An important component of the SWAT model is snowmelt hydrology in watersheds where river runoff is strongly correlated with spring snowmelt (Pradhanang et al. 2011; Abbas et al. 2019), such as in Alpine catchments (Grusson et al. 2015; Tuo et al. 2018). The best simulation of runoff from snowy catchments depends upon how well precipitation processes are modelled within a hydrological model (Troin & Caya 2014). For instance, hydrological models that have been developed to simulate runoff as a function of rainfall (Samadi et al. 2019) are not suitable for snowy catchments. Therefore, an additional snowmelt module is incorporated in the SWAT model (Martinec et al. 2008).

The SWAT model uses the temperature index approach along with elevation bands to simulate streamflow from snow and glacier-fed catchments. In the SWAT, a snowmelt module is incorporated which approximates both snow and glacier melt. SWAT has the capability to model both snow/glacier melt and rainfall–runoff processes but lacks an explicit glacier module. Numerous studies (e.g. Jain 2015; Babur et al. 2016; Garee et al. 2017; Abbas et al. 2019; Anjum et al. 2019) have shown that the SWAT model has the capability to simulate streamflow efficiently in snow/glacier-fed catchments. In this study, snow and glacial melt processes are simulated by temperature index along with elevation bands approach with automatic calibration tool SWAT-CUP. Previous studies revealed that setting up elevation bands and snowmelt parameters, the SWAT model has successfully simulated streamflow in glaciered watersheds in the UIB (Garee et al. 2017; Sarah et al. 2021) and other parts of the world (Rahman et al. 2012; Omani et al. 2017).

The model was set up as follows: (i) delineation of watershed and sub-watershed, (ii) definition of hydrological response units (HRUs) and elevation bands, (iii) model run and parameterization and (iv) calibration, validation and uncertainty analysis.

The delineation of watersheds into sub-basins and HRUs was implemented using ArcSWAT in the ArcMAP10.5 environment. The main river network and tributaries were generated from the SRTM DEM. The UIB was divided into three watersheds, namely the KBR at Nowshera outlet just before the confluence of Kabul River and Indus River, the UIRB at Tarbela outlet and the UJRB at Mangla outlet (Figure 1). SWAT divides the catchment into sub-basins that connect to a river system (Winchell et al. 2007). The three major watersheds of UIRB, UJRB and KRB were further divided into 27, 7 and 11 sub-basins, respectively (Supplementary Material, Figure S8) and then into HRUs, which are the smallest units in watershed delineation using the DEM. For HRU definition, four classes of slope with ranges 0–10, 10–15, 15–20 and above 20% were defined which resulted in 144, 53 and 204 HRUs in the KRB, UJRB and UIRB watersheds, respectively. The water balance components computed at the HRU level are accumulated at each sub-basin level and then carried forward to the outlet of the catchment (Grusson et al. 2015).

Water balance components were simulated for each HRU and then added to each sub-watershed and routed to the outlet of each whole watershed. In the third step, the SWAT model was run based on the prepared CFSR and GCM data. The initial 10-year period (1980–1989) was used as the ‘warm up’ period and streamflow simulated first for the calibration period (1990–2007) and from 2008 to 2013 for validation. In the fourth step, the model was calibrated and validated by using the SWAT-CUP SUFI-2 algorithm (see Section 3.6). The calibrated and validated SWAT model was then run using historical data (1974–2004) and then using data from each GCM output under RCP4.5 and RCP8.5 to simulate future water balance components.

Hydrological models play a vital role in water resource management but various sources of uncertainty make it necessary to implement uncertainty analysis (Zhao et al. 2018). The accuracy of hydrologic model simulation is achieved by parameter sensitivity and uncertainty analysis. Yang et al. (2008) explained various techniques of uncertainties for the better calibration process in the SWAT model in a study conducted in Chaohe basin, China. Calibration and uncertainty analysis are linked together, and calibration results should not be presented without quantification of the degree of uncertainty of model prediction (Abbaspour et al. 2015).

SWAT-CUP is equipped with different tools like Parameter Solution (ParaSol), Generalized Likelihood Uncertainty Estimation (GLUE) and Sequential Uncertainty Fitting (SUFI2) to quantify parameter sensitivity and uncertainty analysis. Among these algorithms, SUFI-2 gained the most popularity which is widely used for parameterization, calibration, validation, sensitivity and uncertainty analysis (Abbaspour et al. 2007; Sao et al. 2020) achieving good prediction uncertainty ranges with fewer runs in large-scale models (Abbaspour et al. 2015). In this study, parameter sensitivity and uncertainty analysis performed by SUFI-2 imbedded into a platform SWAT-CUP (Abbaspour 2012). SUFI-2 is the most widely used method for carrying out calibration, validation, sensitivity and uncertainty analysis (Abbaspour et al. 2007).

For parameter optimization, the SUFI-2 algorithm is used by importing SWAT model simulations into SWAT-CUP (Abbaspour et al. 2007). The SUFI-2 algorithm is used for parallelization with a large number of parameters and observed data from many outlets simultaneously. The SUFI-2 algorithm brackets most of the observed data within a range of 95% prediction uncertainty (95 PPU) band (Supplementary Material, Figures S5–S7) which represents the model prediction (Abbaspour 2015). In SUFI-2, most of the measured data is bracketed with the smallest possible uncertainty band to get optimum values of p-factor (i.e. >0.7) and r-factor (i.e. <1.5) (Abbaspour et al. 2015). A large set of parameter intervals are used in SUFI-2 so that observed data is bracketed by the 95 PPU and then the parameter interval range is narrowed to reduce the uncertainty interval. The goodness of calibration/uncertainty is measured by p-factor and r-factor. p-factor is the percentage of measured data enveloped by 95 PPU, whereas r-factor is the thickness of the 95 PPU envelope divided by the standard deviation of the observed data to evaluate uncertainty. In this exercise, it used to capture most of the observations in the 95 PPU while at the same time having a small envelope. A simulation, which corresponds exactly to measured data, will have a p-factor of 1 and an r-factor of 0. For each iteration, 300–500 simulations are generally performed until a satisfactory p-factor and r-factor is achieved. After achieving a p-factor close to 0.7 and an r-factor close to 1, the ranges of parameters achieved were considered as a calibrated parameter.

Before calibrating the SWAT model, 23 parameters were selected based on the literature review (Khan et al. 2014; Garee et al. 2017). We performed 2,000 simulations in SUFI-2 as recommended in previous studies (Tuo et al. 2016; Zhao et al. 2018) in order to filter out less sensitive parameters to ensure efficient calibration, a global sensitivity analysis performed in SUFI-2. Global sensitivity analysis helps to rank the most sensitive parameters using t-stat and p-value. The larger absolute value of t-stat represents the most sensitive parameter and the p-value measures the significance of the sensitivity of a parameter. It was observed that CN2 was the most sensitive parameter followed by SFTMP, GW-DELAY, REVAPMN, SMTMP and ALPHA-BF, and the remaining parameters were relatively less sensitive for streamflow simulation. The number of parameters reduced from 23 to 12, 12 and 13 parameters for the UJRB, KRB and UIRB, respectively. The selected parameters included three parameters for surface flow (CN2, SURLAG and SLSUBBSN), one for soil water (SOL-AWC), one parameter on evaporation (ESCO), five parameters on groundwater (ALPHA-BF, GW-DELAY, GWQMN, GW-REVAP and REVAPMN) and five parameters on snow with elevation bands (SFTMP, TIMP, SMFMN, SMFMX and SMTMP).

The model was parameterized by performing at least 500 simulations for each iteration and 5–7 iterations for each watershed. For instance, if the simulated base flow was too high, then changes to deep percolation loss (GWQMN), groundwater-revap coefficient (GW-REVAP) and decreasing threshold depth of water in shallow aquifer (REVAPMN) were made to better fit base flow simulations with the observed flow in the UJRB. The Nash–Sutcliff efficiency (NSE) is one of the most commonly used criteria for the calibration and validation of hydrological models and describes the goodness of fit of the observed and simulated plots (Nash & Sutcliffe 1970). Its value ranges from 0 to 1, with values >0.5 considered acceptable and 1 indicating a perfect match between the observed and simulated data. The coefficient of variation (R²) values ranges from 0 to 1, with values higher than 0.5 considered acceptable. PBIAS is an error index, which calculates the average tendency of the simulated data to be either larger or smaller than their measured counterparts, with an optimal value of 0 (Gupta et al. 1999). Both the R² and Nash–Sutcliffe calibration and validation values obtained for the watershed outlets of Tarbela, Mangla and Nowshera were ranked as good to very good (Table 3). The percent bias (PBIAS) is an average tendency of simulation data to be larger or smaller than the observed data (Gupta et al. 1999). Here, the simulation results in the sub-basins of the UIB were good to very good by this performance measure. The magnitude of peak flows in the UJRB and KRB remained the lowest in the calibrated periods during 1999–2002, while lower flows in this particular period were not evident in the UIRB. The underestimation of monthly flows was simulated in some years of the UJRB and KRB while overestimation in the UIRB. However, the monthly timing of peak flows was simulated with good accuracy, which suggests that the SWAT model could be used for further investigation of hydrological components in the studied basins (Figure 2).

The SWAT model was evaluated by using three of the five performance evaluation parameters available in the SUFI-2 algorithm including NSE, R² and PBIAS (Table 4).

Figure 3 and Supplementary Material, Table S7 depict annual and seasonal basin-wide precipitation changes over the three sub-basins (UIRB, UJRB and KRB) of the UIB, respectively. Annual precipitation projections were predicted by using five GCMs (GFDL-ESM2M as GCM1, HadGEM2-ES as GCM2, IPSL-CM5-LR as GCM3, MIROC as GCM4 and NoerESM1-M as GCM5) under two RCPs (4.5 and 8.5) for annual and seasonal (winter [DJF], spring [MAM], summer [JJA] and autumn [SON]) precipitation changes in the near future (2020–2050) and end of the century (2070–2099). Results are presented as average annual and seasonal increases or decreases in precipitation (percentage change) for the near future and end of the century with respect to the reference period (1974–2004).

In the UIRB, GCM1 and GCM2 projected an increase in the average annual precipitation (+6 to +9.1% and +8.1 to +12.2%), while GCM3, GCM4 and GCM5 projected a decrease in the near future and end of the century under RCP4.5 and RCP8.5. However, the average of five GCMs assigned as MMM was projected to increase (+4.5%) under both emission and future scenarios.

Greater variations were found in the seasonal precipitation patterns. GCM1 and GCM2 projected an increase in precipitation (+20 to +35%), whereas oneGCM3 projected a decrease (−18.3%) in precipitation for all four seasons under both emission scenarios and for both periods. The highest increase in precipitation (+16.7%) was projected during summer for the UIRB.

In the UJRB, GCM2 and GCM3 predicted an increase in the average annual precipitation for the reference period (1974–2004), near future, and by the end of the century under both emission scenarios. However, three GCMs (GCM1, GCM4 and GCM5) projected a slight decrease in precipitation with respect to the reference period, but at the end of the century, an increase in monsoon precipitation was projected (+8 to +71%). A consistent increase in winter (DJF) precipitation was observed under RCP4.5 in the near future and end of the century, while a significant decrease in winter precipitation was projected by all GCMs under RCP8.5 in the near future. The average annual precipitation was observed to increase under RCP4.5 in the near future, while no significant change was observed at the end of the century. The MMM projected an increase in monsoon precipitation but a decrease in winter precipitation in the near future under both emission scenarios. Overall, MMM projected an increase in precipitation under both emission and future scenarios over the UJRB. The highest increase in precipitation was projected during summer by MMM under both emission scenarios and future scenarios from the UJRB.

In the KRB, no significant change in annual precipitation was projected in the near future (+6 to −5.2%) and at the end of the century (+12 to −19.5%). Some individual GCMs projected a significant decrease in the average annual precipitation (e.g. GCM3 by −19.5%). Similarly, the average annual precipitation projected by MMM decreased by −7.5%, while a significant increase (+28.5%) in monsoon precipitation was projected under RCP8.5 in the near future and end of the century.

In the sub-basins of the UIB, monsoon precipitation was projected to increase (+7 to +36.8%) by MMM under RCP4.5 and RCP8.5 at the end of the century. In all three sub-basins, GCM precipitation projections varied considerably. However, MMM precipitation projections show a consistent increase in the summer monsoon in the UIRB, UJRB and KRB but a decrease in the average annual precipitation in the KRB. This decrease in the average annual precipitation was due to higher decreases in precipitation during winter, spring and autumn (Supplementary Material, Table S1).

Major properties of the downscaled projections of maximum temperature (T_max) and minimum temperature (T_min) for annual and seasonal projections over the UIRB, UJRB and KRB are presented in Figures 4 and 5.

The downscaled and bias-corrected projections from the selected GCMs showed an increase in T_max ranging from 1.2 to 4.1 °C under RCP4.5 and 1.9 to 7.5 °C under RCP8.5. Similarly, the increase in T_min ranged from 1.3 to 4.4 °C under RCP4.5 and 1.7 to 8.2 °C under RCP8.5.

In the UIRB, an increase in T_max (2.3–6.3 °C) was projected under RCP8.5 in the near future and end of the century (Figure 4). Similarly, the highest increase in T_min was 2–6.3 °C under RCP8.5 in the near future and end of the century. The increase in T_max and T_min under both emission and future scenarios was not much different over the UIRB in the near future and end of the century, with the highest increases in the spring.

On annual basis, the highest average annual increase in T_max and T_min was projected from 2.9 to 3.1 °C, respectively, in the near future in the UJRB. The highest increase in T_max and T_min was projected by GCM4 (8.1 to 8 °C), while the highest projected increase in T_max and T_min by MMM was 5.8 to 5.6 °C, respectively, by the end of the century. On a seasonal basis, the highest rise in T_max and T_min was projected to be 3.9 °C, while 9.5 to 10 °C during spring by the end of the century under both emission scenarios.

In the KRB, the highest rise in average annual T_max and T_min in the near future was projected to be from 2.6 to 3.1 °C, while 4.6 to 8 °C by the end of the century under both emission scenarios by GCM4. On a seasonal basis, the highest rise in T_max and T_min was projected to be between 3.2 °C under RCP4.5 and 3.9 °C under RCP8.5 during autumn in the near future, and projected an increase in T_min and T_max from 9.5 to 10.1 °C under RCP8.5 during the spring by the end of the century. However, the highest rise in T_max (6.1 °C) and T_min (5.8 °C) was projected during winter for both emission and future scenarios. The results suggest higher warming during the winter by the end of the century over the KRB.

Overall, the MMM projected an average annual rise in T_max and T_min in all three sub-basins, while an increase in precipitation in the UIRB under RCP4.5 and RCP8.5 in the near future and end century. The maximum rise in T_max was projected in the KRB, while the highest increase in T_min was projected in the UIRB both in the near future and by the end of the century. The average annual precipitation (Figure 3) was projected to increase in the UJRB in the near future but decrease in precipitation in the KRB under both emission and future scenarios.

The outputs from the five GCMs under two representative concentration pathways (RCP4.5 and RCP8.5) in the near future and end of the century were used to assess the impact on future water balance components by using the SWAT model. Simulated changes in hydrological components with reference to the historical period (1974–2004) in the near future and end of the century in sub-basins (UIRB, UJRB and KRB) were analysed both annually and seasonally.

The average annual and seasonal water yield and evapotranspiration were projected to increase with reference to the historical period both in emission and future scenarios (Supplementary Material, Table S3). In the near future, all GCMs projected peaks in WYLD during July, while by the end of the century, these peaks were more widely spaced from May to July (Figure 6). This shift in water yield indicates an early melting of snow and ice in the year by the end of the century.

All GCMs projected higher increase in water yield during April (41%), May (65%) and June (43%) with reference to the historical period (1974–2004). All GCMs show increase in seasonal water yield but there is a range of values in GCMs outputs for each future period. For instance, the highest increase in water yield was projected by GCM4 (162%) while the lowest increase was projected by GCM2 (9%) during May in the near future under RCP8.5. At the end of the century, GCM1 and GCM4 projected water yield by 15 and 489% during May by GCM1 and GCM4 under RCP4.5 and RCP8.5 respectively. Annual water yield anomalies vary from increases of 2 to 49% by GCM2 and GCM4 respectively in the near future under RCP8.5 and from increases of 23 to 128% by GCM1 and GCM2 under RCP4.5 and RCP8.5 respectively at the end of the Century. Most of the GCMs projected an increasing trend in annual water yield in the near future but a decreasing trend by the end of the century for the UIRB (Figure 7). Under warming in the UIRB, snow and glacier melting will accelerate in the near future, whereas the shrinking glacier area will result in reduced water supplies at the end of the century.

Average annual and seasonal soil moisture was projected to decrease with reference to the historical period (1974–2004) in the UIRB in the near future and end of the century under both emission scenarios. Similarly, mean monthly soil moisture was projected to decrease in the near future and end of the century under both emission scenarios (Supplementary Material, Table S3).

The average annual and seasonal water yield and soil moisture were projected to decrease, while evapotranspiration was projected to increase both in emission and future scenarios (Supplementary Material, Table S4). The highest decrease in water yield was projected during summer. However, the average annual and seasonal evapotranspiration was projected to increase with reference to the historical period both in the near future and end of the century. The mean monthly water yield projected by each GCM over the UJRB showed great variations (Figure 8). Mostly, GCMs projected peaks of water yield twice a year during May and July in the UJRB. The first peak of water yield was due to the melting of snow, while the second peak in water yield was due to monsoon precipitation in the UJRB. An early shift (from May to April) in water yield was projected towards the end of the century. On a seasonal basis, water yield was projected to decrease during spring, summer and autumn, while water yield was projected to increase during winter under both emission and future scenarios (Supplementary Material, Table S4).

GCM1 and GCM5 showed a decreasing trend in water yield under both emission scenarios in the near future and end of the century, respectively (Figure 9). However, GCM2 and GCM3 showed an increasing trend in water yield towards the end of the century. Overall, the average water yield was projected to remain unchanged by the MMM projection.

The average annual and seasonal water yield and soil moisture were projected to decrease both in emission and future scenarios (Supplementary Material, Table S5). The highest increase and decrease in seasonal water yield were projected during winter (24%) and summer (−69%) by the end of the century. This increase in water yield during winter shows a shift from frozen to liquid precipitation due to a rise in temperature over the KRB.

Most of the GCMs projected a decrease in water yield except January, February and March when some of the GCMs projected an increase in water yield (Supplementary Material, Table S4). In the near future, changes in monthly flows are projected to vary from −46% decrease to +51% increase in water yield during January and September by GCM3 and GCM1, respectively, under RCP8.5. At the end of the century, changes in monthly flows are projected to vary from −59% decrease to +65% increase by GCM3 and GCM1, respectively, under RCP8.5. These results indicate that projected changes in water yield are more widespread among GCMs at the end of the century than in the near future. In the near future, annual water yield anomalies vary from −26 to +9% by GCM3 and GCM2 under RCP4.5 and RCP8.5, respectively. At the end of the century, anomalies for water yield vary from −38% decrease to +2% increase by GCM1 and GCM2 under RCP4.5 and RCP8.5, respectively.

Mean monthly water yield showed an early shift in peaks from May in the near future to April by the end of the century (Figure 10). Trends in water yield from the KRB projected by each GCM varied greatly under both emission and future scenarios (Figure 11). For instance, GCM1 and GCM4 projected a significant decreasing trend in water yield in the near future under both emission scenarios. Similarly, GCM4 and GCM5 projected significant decreasing trends in water yield by the end of the century under RCP4.5 and RCP8.5 respectively. However, the average water yield projected by MMM showed no significant trend under both emission and future scenarios.

In the near future, all GCMs projected the highest increase in evapotranspiration during May (Supplementary Material, Table S5). However, some GCMs projected peaks in evapotranspiration during April, while others during May by the end of the century under RCP4.5 and RCP8.5. These results clearly indicate a shift in evapotranspiration from May to April by end of the century, which is likely to affect snowmelt discharge in the KRB.

Future climate change will result in an increase in annual water yield from the nival–glacial regime (i.e. UIRB) in the near future. However, annual water yield is projected to decrease from the pluvial–nival regime (i.e. UJRB and KRB) in the near future as well as by the end of the century. Seasonal water yield is projected to increase during spring and summer from the UIRB under both emission and future scenarios. However, the largest decrease in seasonal water yield is projected during the spring and summer seasons from the pluvial–nival regime, which indicates a decrease in solid precipitation during winter and an increase in evapotranspiration during the summer under the modelled emission scenarios.

Pakistan is highly vulnerable to climate change. For the sustainable use of water resources, it is important to be able to predict climate change impacts on future hydrology and quantify water balance components such as water yield, evapotranspiration and soil moisture at a watershed scale.

The results of this study from precipitation projections in the sub-basins of the UIB showed wide variation which makes the future predictions of water availability in the UIB highly uncertain in the long run (Supplementary Material, Tables S7–S9). For instance, the annual precipitation projected over the UIB ranges from +5 to +25% (Akhtar et al. 2009; Immerzeel et al. 2010; Forsythe et al. 2014; Lutz et al. 2016), while others projected annual precipitation ranging from +40 to −20% (Chaturvedi et al. 2014; Lutz et al. 2014; Khan et al. 2015; Hassanyar & Tsutsumi 2017; Haider 2018; Hayat et al. 2019) over the UIB. Large variations in precipitation projections and trends in the study basins were quite similar to the analysis by Palazzi et al. (2015) from 32 CMIP5 GCMs over the HKH region. These results show that the spread of precipitation projected by each GCM over the UIB was large due to the complex climate and topography of the region. Although subject to large variation, discrepancies in inter-model projections were greater for precipitation than for temperature.

The summer monsoon is a major contributor to river inflow from the pluvial–nival regime, whereas the predicted summer monsoon precipitation showed a remarkable spread over the UIB. The intensity in the monsoon precipitation was projected to increase over the UJRB (Palazzi et al. 2013; Babur et al. 2016; Lutz et al. 2016; Azmat et al. 2018; Bokhari et al. 2018; Sidiqi et al. 2018; Jury et al. 2020). The spread in monsoon precipitation was predicted over the KRB. For instance, GCM2 and GCM3 projected +91% increase to −48% in monsoonal precipitation. Our results of monsoon precipitation are consistent with earlier studies (Palazzi et al. 2013; Rajbhandari et al. 2014; Babur et al. 2016; Bokhari et al. 2018; Sidiqi et al. 2018). The contradictory results of precipitation among GCMs make the water availability highly uncertain in the long run (Lutz et al. 2014). However, it is pertinent to mention here that the MMM-based projections give more confidence in reproducing changes in precipitation than using a single GCM in predicting precipitation over the UIB.

Numerous studies have predicted temperature changes at a greater rate in the UIB than the global average warming (Forsythe et al. 2014; Rajbhandari et al. 2014; Ali et al. 2015; Lutz et al. 2016; Azmat et al. 2018; Jury et al. 2020). Similarly, our analysis projected warming in the UIB which was higher than other parts of the world. The projected average warming by the end of the century ranges from 3.4 to 6.3 °C in the UIRB, 3.3 to 5.7 °C in the UJRB and 3.6 to 6.2 °C in the KRB (Supplementary Material, Table S2). The global average annual warming projected by Knutti & Sedláček (2013) was 1.8 to 4.4 °C under the same RCPs (i.e. RCP4.5 and RCP8.5) by the end of the century. A rise in temperature during the winter and spring is alarming for the hydrological system with implications for snow accumulation and snow melting in this region (Figures 4 and 5). This rise in temperature during winter and spring is consistent with earlier studies (Azmat et al. 2018; Hasson et al. 2019). The increase in temperature during the winter may reduce frozen precipitation leading to a reduction in snow accumulation, which will ultimately lead to reduced snowmelt availability for summer discharge. Similarly, the rise in temperature during the spring may trigger early melting of snow leading to an early shift in peak flows from the sub-basins of the UIB. However, there is no agreement on results regarding the sensitivity of runoff to temperature in previous studies (Jung et al. 2013; Jepsen et al. 2016; Wang et al. 2016; Zhang et al. 2020). In the UIB, future snowmelt is expected to decrease slightly due to increasing temperature and limited change in precipitation projected by all ensemble members in both RCPs, while a shift towards more liquid precipitation due to warming in winter (Lutz et al. 2014). According to Bokhari et al. (2018), annual and seasonal warming in the UIB may negatively affect snow accumulation during the winter and has the potential to accelerate glacier melting during the summer. The projected changes in river flows from the Kabul and Indus Rivers reveal the difference in responses to climate change in the sub-basins of the UIB (Lutz et al. 2014). Water yield was projected to increase in the near future and by the end of the century from the UIRB, while water yield was projected to decrease from the KRB (Figures 6 and 10).

It is well documented in previous studies that rising temperatures generally increase actual evapotranspiration, which results in a decrease in water yield and soil moisture (Sun et al. 2016; Guo et al. 2017; Liu et al. 2019). According to Singh & Bengtsson (2005), snowy catchments are more sensitive to the effect of warming in the reduction of water yield due to an increase in evaporation and a decrease in snow accumulation. However, for complex catchments such as UIB, a decrease in melt from seasonal snow can be counterbalanced by an increase in melt from glaciers. Therefore, it could be inferred from the complex hydrology of the UIRB dominated by snow and glaciers that an increase in water yield was a result of higher warming under both emission and future scenarios (Table 4). Basins with high snowmelt contributions (UJRB and KRB) showed a reduction in water yield under warming conditions due to increases in evapotranspiration (ET) for future scenarios (Supplementary Material, Tables S4 and S5).

Future projections of water yield at three outlets of the sub-basins of the UIB showed highly contrasting hydrological signals. A significant decrease in annual water yield from the UJRB and KRB could be due to higher projected evapotranspiration (Supplementary Material, Tables S4 and S5). Similarly, Lutz et al. (2016) and Dahri et al. (2021) projected a decrease in river inflows at two gauging stations (e.g. UJRB at Mangla and KRB at Nowshera) but an increase in water yield from the UIRB (Hasson et al. 2019; Khan et al. 2020; Dahri et al. 2021). Our results showed a consistent increase in seasonal inflows from all three sub-basins during spring which is in line with the results of Dahri et al. (2021) in the KRB. The increase in water yield was only predicted during winter which indicates that future winter warming may result in an increase in precipitation in liquid form in the KRB (Supplementary Material, Table S5). However, a decrease in water yield during summer months was consistent with others studies in the KRB (Babur et al. 2016; Lutz et al. 2016; Iqbal et al. 2018).

According to Lutz et al. (2014), the projected water availability from the UIB was very uncertain with a range from −15 to +60%, while Hasson et al. (2016) projected a 220% increase by the end of the century. Immerzeel et al. (2010) projected a decrease in upstream water supply (−8.4%) for the UIB. In the lower altitude sub-basins of the UIB (UJRB and KRB), rainfall–runoff and snowmelt components dominate river discharge from the sub-basins of the UIB. The projected flows from the UIB showed by some models to decrease as much as −50% due to decreases in glacial melt, while others predicted an increase in flow due to higher rates of snowmelt, precipitation and temperature (Ragettli et al. 2013). Despite large differences between basins (i.e. UIB, Ganges, Brahmaputra and Salween) in the Himalayas and between tributaries within basins, Lutz et al. (2016) projected an increase in water availability from the UIB primarily due to accelerated glacier melt until 2050. Nie et al. (2021) projected an increase in streamflow by 51% by the end of the century under RCP8.5, while a 40.3% increase in water yield was projected from the UIRB in this study (Supplementary Material, Table S3).

In the UIRB, a consistent increase in water yield is projected by all GCMs in the near future, while a decline in peak flows towards the end of the century (2070–2100). These results show projected peaks in water yield around 2070 under RCP8.5 (Figure 7), while Huss & Hock (2018) also projected peak water from the UIB around 2070. This trend in peak flows is more consistent under RCP8.5 than RCP4.5. These results show an agreement with Nie et al. (2021), where they predicted a delay in glacier melting and a decline in peak flows late in the 21st century under RCP8.5 (2,064 ± 19). This decline in peak flows may be attributed to declining glaciered runoff towards the end of the century from the UIRB. Our analysis contrasts with Hasson et al. (2019) who predicted a 60% increase in annual and unchanged peaks in the KRB by the end of the century.

Monthly peak flows from the UIRB were projected to shift to early in June from July towards the end of the century (Figure 6). This shift in water yield in the highly glaciered catchment of the UIRB showed an early onset of snow and glacier melt. This was most likely due to an increase in temperature during spring. Seasonal streamflow is predicted to rise in June by the end of the century, indicating a shift in peak flows in the UIB (Huss & Hock 2018; Nie et al. 2021). An early start of snow melting from the glaciered sub-basins of the UIB will lead to exposed glacier ice surfaces which themselves will then start melting earlier in the year. Monthly peak flows were projected to increase significantly during April–July months with an overall annual increase of 47–58% in the near future and end of the century under RCP4.5 and RCP8.5, respectively. Dahri et al. (2021) also projected strong increases in Indus-Tarbela inflows (i.e. 17–73.6%) under MIROC5 ensembles. Lutz et al. (2016) anticipated no change at downstream Tarbela except slight increases for the upstream Karakoram river inflows. Hasson et al. (2019) also projected increases in downstream Tarblela inflows, while no shift in peak flows which is likely to shift from June to July under RCP8.5 at the end of the century (Figure 6). A similar shift in WYLD was detected in the UJRB where peaks in water yield were shifted towards April from May at the end of the century. These results are consistent with other studies over the UJRB (Babur et al. 2016; Mahmood & Jia 2016; Azmat et al. 2018). Small differences in the projection of WYLD could be due to differences in selecting GCMs, hydrological models, historical periods and future scenarios. Similarly, shifts in the timing and magnitude of snowmelt runoff have been studied across the HKH region (Tahir et al. 2011; Lutz et al. 2016; Khan et al. 2020). A shift in the timing of streamflow hydrograph by almost 6 weeks could be a potential impact of future atmospheric warming depending on changes in snowfall–rainfall segmentation and snowmelt timing (Costa-cabral et al. 2012). Our results agreed well with those of previous studies (e.g. Lutz et al. 2016; Huss & Hock 2018; Dahri et al. 2021; Nie et al. 2021) such as peak water around 2070 and rise in WYLD in June from the UIRB, which is highly glaciered. However, changes in water yield from the pluvial–nival regime (i.e. UJRB and KRB) predict otherwise. However, predicted water yield from the pluvial–nival regime (i.e. UJRB and KRB) shows a decreasing trend towards the end of the century.

This study assesses future climate change projections over the UIB (i.e. UIRB, UJRB and KRB) and the consequences of these changes for water resources. Temperature and precipitation datasets from CMIP5 GCMs were used to force the SWAT model for the simulation of hydrological regimes in the sub-basins of the UIB under two emission and future scenarios.

The main findings revealed a considerable increase in annual and spring temperature with varying magnitudes in precipitation. The projected changes in temperature and precipitation from climate models vary widely at the sub-basin scale. A higher spread in precipitation projections makes the future hydrological predictions very uncertain over the UIB. Due to the complex climate system of the region, accurate prediction is still difficult due to poor capture of this complexity in models. These variations are further reflected in the simulated flows from sub-basins of the UIB. For instance, the nival–glacial regime is temperature driven due to the dominant component of snow and glacier melting. Therefore, variations in the magnitude of projected changes in water balance components from five GCMs and scenarios are smaller. On the other hand, the projected changes in water balance from five GCMs and scenarios have shown large variations from the precipitation-driven river flow in sub-basins of the UIB (e.g. UJRB and KRB). A large spread in precipitation results was noted caused by variances in predicted precipitation among GCMs. However, variation is addressed considerably by using MMM in the simulation of water balance for future scenarios. Therefore, the use of a single GCM for hydrological simulations for future scenarios is not recommended over the whole UIB.

The implications of future warming on snow/glacier melting are clear from the simulations of water balance components across the UIB. Warming during winter is likely to increase the proportion of rainfall compared to snowfall and thus reduce the snowpack, while early melting of snow during spring will greatly affect the timing of peak flows (earlier and less) from the sub-basins of the UIB. Similarly, a significant increase in water yield in the UIRB may result from higher rates of snow/glacier melting within highly glaciered sub-basins of the UIB. Future research should look towards incorporating variation in ice/glacier extent in the region explicitly into hydrological models, which at present is a current limitation of the SWAT model used here.

The predicted changes in peak flows from the UIB have implications for hydroelectricity, industry, downstream irrigation and crop patterns, as well as society in general. Therefore, a better understanding of climate change impacts on water resources is essential for decision-making in water-dependent sectors. In the future, policy-makers and stakeholders should keep in mind the potential upcoming changes in streamflow such as early snow melting in the KRB and UJRB and higher stream flows from the UIRB. The outcomes of this research are crucial for both scientific community and policy-makers in Pakistan and within the wider region for understanding climate change impacts on water resources of the UIB. In this regard, GCM-based predictions can be utilized to bridge the knowledge gap for near- and longer-term changes in the hydrology of the UIB.

The authors are grateful to the WAPDA for providing data on river flow at the Mangla outlet of the UJRB. This work was completed as part of PhD thesis which was funded by the Punjab Higher Education Commission (PHEC).

The authors declare that there is no conflict of interest.

Data cannot be made publicly available; readers should contact the corresponding author for details.