Impacts of Climate Change on the Water Availability, Seasonality and Extremes in the Upper Indus Basin (UIB)

Projecting future hydrology for the mountainous, highly glaciated upper Indus basin (UIB) is a challenging task because of uncertainties in future climate projections and issues with the coverage and quality of available reference climatic data and hydrological modelling approaches. This study attempts to address these issues by utilizing the semi-distributed hydrological model “Soil and water assessment tool” (SWAT) with new climate datasets and better spatial and altitudinal representation as well as a wider range of future climate forcing models (general circulation model/regional climate model combinations (GCMs_RCMs) from the “Coordinated Regional Climate Downscaling Experiment-South Asia (CORDEX-SA) project to assess different aspects of future hydrology (mean flows, extremes and seasonal changes). Contour maps for the mean annual flow and actual evapotranspiration as a function of the downscaled projected mean annual precipitation and temperatures are produced and can serve as a “hands-on” forecast tool of future hydrology. The overall results of these future SWAT hydrological projections indicate similar trends of changes in magnitudes, seasonal patterns and extremes of the UIB—stream flows for almost all climate scenarios/models/periods—combinations analyzed. In particular, all but one GCM_RCM model—the one predicting a very high future temperature rise—indicated mean annual flow increases throughout the 21st century, wherefore, interestingly, these are stronger for the middle years (2041–2070) than at its end (2071–2100). The seasonal shifts as well as the extremes follow also similar trends for all climate scenario/model/period combinations, e.g., an earlier future arrival (in May–June instead of July–August) of high flows and increased spring and winter flows, with upper flow extremes (peaks) projected to drastically increase by 50 to >100%, and with significantly decreased annual recurrence intervals, i.e., a tremendously increased future flood hazard for the UIB. The future low flows projections also show more extreme values, with lower-than-nowadays-experienced minimal flows occurring more frequently and with much longer annual total duration.


Introduction
Large hydrological changes are expected in the coming decades throughout the world, mainly due to the changing climate with its unequivocal rise in global temperatures [1,2]. The latter will result in an accelerated water cycle, accompanied by an increasing probability of the occurrence of extreme hydrometeorological events and, possibly, altered seasonal patterns in different regions of the world [3,4]. The hydrological impacts of climate change may vary considerably in magnitude, 1. to minimize uncertainties in climate change hydrological impact study in the UIB by utilizing a set of improved reference climate data for the setup and calibration of the hydrological model; 2.
to simulate and project future flow regimes in UIB under five different climate change scenarios using the SWAT hydrological model; and 3.
to assess the impacts of climate change on the water availability, seasonality and extremes in the Upper Indus Basin (UIB)

Study Area-The Upper Indus River Basin (UIB)
The Indus River, one of the largest rivers in Asia with a total length of about 2880 km, has a total drainage area of about 912,000 km 2 , which extends across portions of India, China, Pakistan and Afghanistan [13]. The section of the Indus upstream of the Tarbela Dam, Pakistan (Figure 1), Sustainability 2020, 12, 1283 4 of 27 which, according to our findings, has a length of about 1150 km and drains an area called the UIB, is approximately 165,400 km 2 .
Being a high mountain region, the UIB contains the largest area of perennial glacial ice cover of approximately 15,062 km 2 outside the polar regions, with a volume of around 2174 km 3 of the total ice reserves [35]. Some estimations [36] show even greater glacier cover, attaining an area of up to 12% of UIB (above 19,000 km 2 ). The altitude of the UIB ranges from as low as 455 m to a maximum of 8611 m, and as a result, the climate varies greatly within the basin [10].
The summer monsoon has little effect on the basin, as almost 90% of its area is located in the rain shadow of the Himalayan belt [24]. Hence, except for the south-facing foothills, the intrusion of the Indian-ocean monsoon is limited by the mountains, so that its influence weakens northwestward [37]. Subsequently, the climatic controls in the UIB are quite different from that in the Himalayas on the eastern side. In fact, over the extent of the UIB, most of the annual precipitation originates from the mid-latitude western disturbances and mostly in solid form during winter and spring [25,[38][39][40]. Occasional rains are brought in also by the monsoonal incursions into trans-Himalayan areas [38,40], but even during the summer months, the trans-Himalayan areas do not receive all precipitation from monsoon sources [41].
Sustainability 2020, 12, x FOR PEER REVIEW 4 of 30 according to our findings, has a length of about 1150 km and drains an area called the UIB, is approximately 165,400 km 2 . Being a high mountain region, the UIB contains the largest area of perennial glacial ice cover of approximately 15,062 km 2 outside the polar regions, with a volume of around 2174 km 3 of the total ice reserves [35]. Some estimations [36] show even greater glacier cover, attaining an area of up to 12% of UIB (above 19,000 km 2 ). The altitude of the UIB ranges from as low as 455 m to a maximum of 8611 m, and as a result, the climate varies greatly within the basin [10].
The summer monsoon has little effect on the basin, as almost 90% of its area is located in the rain shadow of the Himalayan belt [24]. Hence, except for the south-facing foothills, the intrusion of the Indian-ocean monsoon is limited by the mountains, so that its influence weakens northwestward [37]. Subsequently, the climatic controls in the UIB are quite different from that in the Himalayas on the eastern side. In fact, over the extent of the UIB, most of the annual precipitation originates from the mid-latitude western disturbances and mostly in solid form during winter and spring [25,[38][39][40]. Occasional rains are brought in also by the monsoonal incursions into trans-Himalayan areas [38,40], but even during the summer months, the trans-Himalayan areas do not receive all precipitation from monsoon sources [41]. Climate variables are usually strongly influenced by topographic altitude. Thus, the northern valley floors of the UIB are arid, with an annual precipitation of only 100-200 mm which increases to 600 mm at a 4400-m altitude, and for higher elevations, glaciological studies suggest annual accumulation rates of 1500 to 2000 mm at 5500 m [40]. The average annual snow cover area in the UIB ranges from 10% to 70%, with its minimum in the summer snow melt period (June-September) and its maximum in the winter snow accumulation period (December-February) [10].
Stream flow in the UIB is generated by the combination of storm runoff in its lower part and snow and glacier runoff in its higher parts [12,13] Nearly half of the total annual flow in the whole Indus basin is contributed by the UIB, with a 86-88% contribution during the summer and only around 12-14% during the winter season [11,32]. Climate variables are usually strongly influenced by topographic altitude. Thus, the northern valley floors of the UIB are arid, with an annual precipitation of only 100-200 mm which increases to 600 mm at a 4400-m altitude, and for higher elevations, glaciological studies suggest annual accumulation rates of 1500 to 2000 mm at 5500 m [40]. The average annual snow cover area in the UIB ranges from 10% to 70%, with its minimum in the summer snow melt period (June-September) and its maximum in the winter snow accumulation period (December-February) [10].

Data Description
Stream flow in the UIB is generated by the combination of storm runoff in its lower part and snow and glacier runoff in its higher parts [12,13] Nearly half of the total annual flow in the whole Indus basin is contributed by the UIB, with a 86-88% contribution during the summer and only around 12-14% during the winter season [11,32]. The in situ meteorological observations in the UIB are sparse and mostly taken at (lower altitude) valley stations with very low spatial coverage. In addition, very little data is available for high altitudes. Furthermore, the orography of the UIB region is complex, leading to coactions of different hydro-climatic regimes which affect the amounts, spatial patterns and seasonality of the precipitation.
To cater for the non-representativeness and reported underestimations of the precipitation amounts by all the available data products [29,42], the current study has utilized a new corrected climate dataset generated by Khan and Koch [29,34].
This reference precipitation data set provided by Khan and Koch [29,34] is a combination of three sets of data: (1) 0.125 • by 0.125 • gridded precipitation time series for the 1997-2008 period, generated through the "informed-regionalization" methodology [29] which in turn is based on a relatively dense observing station network; (2) reconstructed precipitation time series for the time period of 1961-1995 for the parts of the UIB inside Pakistan's borders (UIB-Pak) [34]; and (3) modified and corrected APHRODITE daily precipitation data available for the period 1961-2007 for parts of the UIB outside Pakistan's borders [34].
The improvement and correction of the precipitation data has been done through a "Correction and Informed Regionalization (CIR)" approach. In the CIR methodology, the precipitation data is first quality checked and corrected for losses. Then, it is improved for a better spatial coverage through geospatial interpolation, followed by incorporation of a subbasin-specific "precipitation laps rate", based on reported hydro-climatological data and the glacier mass balance dynamics in the UIB.
The temperature data used is acquired from Khan and Koch [34]. Similar to the precipitation dataset, the temperature data (1996 onwards) is also a combination of the spatially interpolated version of the available observed station data with adjustments by the region-specific, elevation-dependent lapse rate. Using the observed monthly time series variations, the temperature data was extrapolated to reconstruct missing data for the period 1961-1996 for the parts of UIB inside Pakistan's border (UIB-Pak) and a modified and corrected version of APHRODITE data for the rest of the UIB outside Pakistan's border. Figure 2 illustrates maps of precipitation and temperatures for the UIB. Full details of the methodology can be explored in References [29,34]. The in situ meteorological observations in the UIB are sparse and mostly taken at (lower altitude) valley stations with very low spatial coverage. In addition, very little data is available for high altitudes. Furthermore, the orography of the UIB region is complex, leading to coactions of different hydro-climatic regimes which affect the amounts, spatial patterns and seasonality of the precipitation.
To cater for the non-representativeness and reported underestimations of the precipitation amounts by all the available data products [29,42], the current study has utilized a new corrected climate dataset generated by Khan and Koch [29,34].
This reference precipitation data set provided by Khan and Koch [29,34] is a combination of three sets of data: (1) 0.125° by 0.125° gridded precipitation time series for the 1997-2008 period, generated through the "informed-regionalization" methodology [29] which in turn is based on a relatively dense observing station network; (2) reconstructed precipitation time series for the time period of 1961-1995 for the parts of the UIB inside Pakistan's borders (UIB-Pak) [34]; and (3) modified and corrected APHRODITE daily precipitation data available for the period 1961-2007 for parts of the UIB outside Pakistan's borders [34].
The improvement and correction of the precipitation data has been done through a "Correction and Informed Regionalization (CIR)" approach. In the CIR methodology, the precipitation data is first quality checked and corrected for losses. Then, it is improved for a better spatial coverage through geospatial interpolation, followed by incorporation of a subbasin-specific "precipitation laps rate", based on reported hydro-climatological data and the glacier mass balance dynamics in the UIB.
The temperature data used is acquired from Khan and Koch [34]. Similar to the precipitation dataset, the temperature data (1996 onwards) is also a combination of the spatially interpolated version of the available observed station data with adjustments by the region-specific, elevationdependent lapse rate. Using the observed monthly time series variations, the temperature data was extrapolated to reconstruct missing data for the period 1961-1996 for the parts of UIB inside Pakistan's border (UIB-Pak) and a modified and corrected version of APHRODITE data for the rest of the UIB outside Pakistan's border. Figure 2 illustrates maps of precipitation and temperatures for the UIB. Full details of the methodology can be explored in References [29,34].  Daily river discharge data for a total of nine (9) hydrometric stations in the study area were collected from the Water & Power Development Authority, Lahore, Pakistan (WAPDA). These hydrometric stations are located on the main river (Indus) as well as its major tributaries (Table 1). It is important to note that, as the population and human intervention in the UIB is minimal, the flow regime is almost undisturbed and natural. As a matter of fact, central and eastern UIB, in particular, has in the past hardly seen any human-induced major influence on the rivers' discharges; meanwhile, there are only minor single-season streamflow diversions for local irrigation, with minor effects in the western UIB [43,44]. There are, however, chances of interventions in the future as at least three major water reservoirs are proposed in the UIB, namely Bhasha dam, Dasu dam and Bunji dam. [45], which would have the capability to influence the flow regime in UIB. The current study has opted for future hydrological simulations without assessing the influences of the proposed dams as the exact timeline of their construction, inauguration or operation is not known. The future climate forcing employed in the current study were adopted from Khan and Koch [17]. The data consisted of bias-corrected versions of five Coordinated Regional Climate Downscaling Experiment (CORDEX)-SA RCM experiments, namely IPSL-CM5A-MR_RCA4, MPI-ESM-LR_RCA4, NorESM1-M_RCA4, CanESM2_RegCM4-4 and iGFDL-ESM2M_RCA4. The downscaling for these five GCMs is done by CORDEX-SA, using two different RCMs (RCA4 and RegCM4). Their RCM outputs are at a considerably finer scale (0.44 • ) than those of the source GCMs.
This data was adopted as an effort to compensate for the uncertainties and variations between the different CMIP5 coupled climate models' (GCM) results (or the regional climate model projections). Khan and Koch [17] recommended using a wider range of future climate scenarios for two representative concentration pathways (RCPs), namely RCP 4.5 and RCP 8.5, to force the SWAT hydrological model. Khan and Koch [17] have carried the selection of possible future climate scenarios through a stepwise procedure, whereby future projections for RCP 4.5 and RCP 8.5 of the r1p1i1 ensemble from all climate models were screened based on three criteria: (1) the range of projected change in the mean, (2) range of projected change in the extremes and (3) skill in reproducing the past climate. The aim of this procedure was then to attain sets of four possible future extreme scenarios groups (wet-warm, wet-cold, dry-warm and dry-cold) for each RCP as well as of one scenario in each RCP group, representing the mean future climate change over the UIB. Using this category approach made it possible to arrive at a limited number of climate models. The final selection of a model was then carried out by assigning ranks based on the weighted score for each of the mentioned selection criteria. Khan and Koch [17] decided to utilize any dynamically downscaled version for the selected GCMs, Sustainability 2020, 12, 1283 7 of 27 ranked at the 1st or 2nd positions based on availability. As RCM projections were already available for the selected GCMs, five CORDEX-SA experiments (Table 2), including IPSL-CM5A-MR_RCA4, MPI-ESM-LR_RCA4, NorESM1-M_RCA4, Can ESM2_RegCM4-4 and GFDL-ESM2M_RCA4, were selected by Reference [17] for bias correction and further processing.
In the CORDEX-SA experiments, the downscaling for these five GCMs is done by using two different RCMs (RCA4 and RegCM4). Their RCM outputs are at a considerably finer scale (0.44 • ) than those of the source GCMs [17].
The five (5) selected (CORDEX-SA) RCM outputs have been further bias corrected using the "distribution mapping technique" [46] for RCP 4.5 and RCP 8.5 for two sets of durations, i.e., mid-century (2041-2070) and end-century (2071-2100). The five (5) selected (CORDEX-SA) RCM outputs for RCP 4.5 and RCP 8.5 are bias corrected for two future periods, i.e., mid-century (2041-2070) and end-century (2071-2100) using the "distribution mapping technique". Major properties of the downscaled projections are given in Table 3. Detailed descriptions of the downscaled and bias-corrected data can be found in Khan and Koch [17].
The downscaled projections by Khan and Koch [17] show changes/increases over all scenarios of temperature, ranging from 2.3 • C to 6.33 • C for RCP 4.5 and from 2.92 • C to as high as 9.0 • C for RCP 8.5. The higher temperature increase in this region may be attributed to the "elevation-dependent warming" phenomenon, whereby a larger temperature increase is expected in the mountains than in the neighboring plains and valleys of the UIB [16,21,22].
The changes in the downscaled and bias-corrected precipitation over the 21st century ranges from a minor increase of 2.2% for the drier scenarios to as high as 15.9% for the wet scenarios. Thus, both temperature and precipitation show increases, as do the extremes, since the probabilities of the wet days are projected to decrease while the precipitation intensities are projected to increase unanimously for both RCPs.
The spatial distribution of the projected changes in temperature and precipitation across the UIB also exhibits certain distinct trends ( Figure 3). Thus, the late-century years (2071-2100) projected that precipitation changes show a decrease in the southeastern parts and an increase in the northeastern parts of the basin, and this holds for both RCPs. The spatial distributions of the projected changes in temperature show similarities for the two RCPs but with a higher increase in the northern and northwestern parts than in the eastern and southern part of the basin.

SWAT Model Description
The SWAT model is a continuous time (long-term yield) process-based semi-distributed hydrological model developed for the US Department of Agriculture (USDA), Agricultural Research Service (ARS), by Dr. Jeff Arnold [49]. It is capable of simulating long-term impacts of any management, climate or vegetation scenarios on the hydrological processes, pollution transport and sediment loading in river basins/watersheds [49,50]. In recent years, SWAT has been used more frequently, especially for the analysis of climate change hydrological impacts on stream flow and water yield [51][52][53].
To model the various physical/hydrological processes mentioned, SWAT requires specific information pertaining to the watershed such as weather, topography, soil properties, land cover, land use and management practices.
In the SWAT model, a river basin or watershed is partitioned into larger subunits called subbasins draining into the stream network and the river system. These subbasins are further divided into a series of smaller units, the hydrological response units (HRUs), which are nonspatial uniform units, each representing unique combinations of soil, land-use and slope. The calculations and simulations of hydrological components, sediment yield and nutrient cycles are first carried for each HRU and then aggregated for the subbasins and routed through the watershed.
The hydrological cycle simulated in SWAT is based on the water balance equation: where SW 0 and SW t are respectively the initial and final soil water content (mm H 2 O) on day i, t is time (days), R day is the precipitation amount reaching the soil surface on day i (mm H 2 O), Q sur f is the surface runoff amount on day i (mm H 2 O), E a is the evapotranspiration on day i (mm H 2 O), w seep represents the interflow on day i (mm H 2 O) and Q gw is the return flow or base flow to the channel on day i (mm H 2 O). The simulated hydrological components in SWAT include infiltration, redistribution, evapotranspiration (ET), lateral subsurface flow, groundwater or return flow, surface runoff, ponds inflow and outflow, transmission losses and water yield [49]. More detailed background and the theories of the SWAT-modelled hydrological process can be found in Reference [54].

Spatial Input Data and SWAT Model Setup
SWAT requires specific information pertaining to the watershed, such as topography, soil properties and land cover/land use, in addition to climate data and management practices in the basin.
The topographic data was provided in the form of a hydrologically conditioned digital elevation model (DEM), acquired from the "Hydrological data and maps based on SHuttle Elevation Derivatives at multiple Scales" (HydroSHEDS) website. This DEM has a resolution of 3 arc-seconds, is void-filled and is a hydrologically conditioned product [55], derived from the Shuttle Radar Topography Mission (SRTM) DEM [56].
The soil data were obtained from the FAO-UNESCO global soil map version 3.6 [57], which provides data for 5000 soil types comprising two layers (0-30 cm and 30-100 cm depth) at a spatial resolution of 5 km. For land use, the "Global Land Cover Characterization" (GLCC) at 1 km spatial resolution [58] was used.
The high-resolution DEM used in the study needed no preprocessing as it was already void-filled and hydrologically conditioned [55] and therefore provided almost no inconsistencies during the stream network delineation. For the watershed delineation process, "ArcSWAT-2012", which is an ArcGIS-ArcView extension and graphical user input interface for the SWAT model, was used, choosing a threshold drainage area of 50,000 ha for the various subbasins. The whole watershed outlet as well as the monitoring gauge stations were manually added, while subbasin outlets were created by the software automatically based on the given threshold. As a result, the study area was finally configured with 173 subbasins, which were further divided into 2825 discrete HRUs based on the different soil and land-use classes, five (5) slope classes and ten (10) elevation bands. The total simulated watershed area in the current model is 165,611 km 2 .

Model Calibration and Validation Setup
The SWAT model was calibrated and validated against daily discharge data individually for each of its five (5) major tributaries (Hunza, Gilgit, Astor, Shigar and Shyok rivers), for parts of UIB (except the tributaries) inside Pakistan's boundary and for UIB (situated in India China and Nepal) covering area upstream of Kharmong gauge station.
The Sequential Uncertainty Fitting SUFI-2 algorithm of the SWAT Calibration and Uncertainty Programs (SWAT-CUP) [59] was employed for parameter optimization. This algorithm is capable of mapping all uncertainties (parameter, inputs, conceptual model, etc.) in terms of parameter ranges as the procedure attempts to cover most of the measured data within the 95% prediction uncertainty (95PPU), which is calculated at the 2.5% and 97.5% levels of the cumulative distribution of all simulated output values. During this process, the user first assigns tangible ranges to a set of calibration parameters, where both (the ranges and the selection of calibration parameter) are guided by literature information, specific knowledge of the study area and the parameters' sensitivity analysis. Once this is done, sets of samples (as many as intended simulations) are drawn from the parameter ranges through Latin hypercube sampling, followed by the SWAT model simulation using each of the sets and processed to evaluate various objective functions, as discussed below.
To quantify the goodness of model performance, i.e., the objective function, for the selected ranges of parameters in terms of calibration/uncertainty levels, two indices, the P-factor and the R-factor, are used. The P-factor denotes the percentage of data that is bracketed by the 95PPU band (ranging from 0 to 1, where 1 shows that all the predictions are within the 95PPU Band), while the R-factor is the average width of the 95PPU band divided by the standard deviation of the measured variable (0 to ∞, with 0 indicating a perfect match) [60][61][62].
For the evaluation of the calibration/validation results, the SUFI-2 algorithm allows the user to select from a range of different objective functions, such as the coefficient of determination (R 2 ), percent bias (PBIAS), Nash-Sutcliffe efficiency (NS) or Kling-Gupta efficiency (KGE). The objective function can easily be reassigned in the post-processing step if required [59]. The current study used NS as the main objective function for the calibration/validation, but the results were also evaluated based on R 2 , PBIAS and KGE as well as the P-factor and the R-factor. Further details of the SWAT model calibration and validation setup are given in Appendices A-F (see Supplementary Materials, Appendices).

Assessment of Changes in Future Hydrology
The calibrated SWAT model was forced with the downscaled future climate projections for two RCPs (RCP 4.5 and RCP 8.5) of the five (5) selected CORDEX-SA climate models (see Section 2.2.3). For each scenario, the future hydrology is assessed for general flow characteristics, seasonality, high and low flows as well as extremes (return periods of extreme discharges) for two future periods, namely mid-century (2041-2070) and late-century (2071-2100).
These future hydrological assessments include calculations of percentage relative changes from the baseline period  for each of the hydrometeorological variables computed. These rates of change are calculated as follows: %Relative change = X p − X r X r × 100 where Xp and Xr denote the means of the hydrometeorological variable over the future projected and the reference periods, respectively. The future changes in glacier/snowmelt contribution are estimated using the following modified hydrological balance equation proposed by Mark and Seltzer [63]: According to which the total volume of water discharging from a catchment (per specified time period), Q t , is equal to the volume of water entering the catchment (Precipitation P) minus losses (i.e., evapotranspiration, E, and aquifer/groundwater recharge, Gw) plus changes in storage ∆g in the system which comprise loss or gain of glacier ice volume, with snow/glacier accumulation (−) and snow/glacier melt contributions (+).
As the losses Gw may be minimal in comparison to the total discharge, especially for large mountainous catchments at an inter-annual time step, this term may be omitted in the equation above. Furthermore, by substituting ∆g with the changes in snow/glacier accumulation/melt contributions (+∆Gm or −∆Gm), one gets the following: In order to assess the future changes in the snow/glacier melt (∆Gm) contribution relative to that in the reference period, the above equation is rearranged as follows: Similarly, future changes in the hydrological extremes are assessed in terms of the changes in the 1st and 99th percentiles of the flows, changes in the low and high flows, along with an analysis of the changes in recurrence intervals of extreme discharges, with return periods of 10, 15, and 30 years. These are calculated from a fitted Log-Pearson III value distribution to the simulated flow maxima for the reference and the two future 30-year periods: 1976-2005, 2041-2070 and 2071-2100, respectively.
In addition to the above analysis, a 2nd order Response Surface Regression model (RSM) is applied to capture the interrelationship of mean annual "Flow" or "ET' with mean annual precipitation and temperature for the projected future hydro-climatology. RSM models are multivariate polynomial models which utilize a set of mathematical equations to designate relationship between the response (e.g., flow) and the independent variables, e.g., temperature and precipitation [64]. Along with the fitted model, contour plots for the response variables "Flow" or "ET' are also generated. The RSM, in addition to giving an overview of flow/evapotranspiration responses for different temperature and precipitation scenarios, is also usable for preliminary predictions of these hydrometeorological variables under diverse future precipitation and temperature changes. The calculations and plot generation have been done utilizing the "Response Surface Analysis" (RSA) program of the NCSS software (NCSS) (e.g., References [65,66]). A brief description of RSM is provided in Appendix H (see Supplementary Materials, Appendices).

Performance of the SWAT Model (Calibration-Validation)
For almost all tributaries included in the SWAT's calibration process, the most sensitive calibration parameters turn out to be the ones related to snowmelt and groundwater and they are usually homogeneous across all the calibrated basins, i.e., there are only minor differences in the final best fitted ranges and values. A list of the calibrated parameters for the various tributary basins is provided in Appendix F (see Supplementary Materials, Appendices).
Results of the calibration and validation of the flow rates for the basin outlet at gauge station Bisham Qila are shown in Figure 4. One can notice that the SWA model is able not only to simulate most of the observed high and low flows at daily and monthly scales but also to capture the seasonal variation of the discharge as well. In fact, as listed in the figure and detailed in terms of the various performance statistical indicators (R 2 , NS coefficient, PBIAS and KGE) for all gauged stations at the tributaries outlets in Table 4, the calibrated SWAT simulations result in high NS, ranging from 0.69 to up to 0.86 for the flows at the tributaries' outlets and between 0.70 and 0.89 at other monitoring points despite the fact that the UIB is a complex watershed to model, particularly on the daily time scale, as it covers a very vast area with a huge glaciated part and diverse hydro-climatic regimes as well as very sparse observation data. Table 4 indicates furthermore that, similar to the calibration period, the statistical results are also very good for the validation period. Results of the calibration and validation of the flow rates for the basin outlet at gauge station Bisham Qila are shown in Figure 4. One can notice that the SWA model is able not only to simulate most of the observed high and low flows at daily and monthly scales but also to capture the seasonal variation of the discharge as well. In fact, as listed in the figure and detailed in terms of the various performance statistical indicators (R 2 , NS coefficient, PBIAS and KGE) for all gauged stations at the tributaries outlets in Table 4, the calibrated SWAT simulations result in high NS, ranging from 0.69 to up to 0.86 for the flows at the tributaries' outlets and between 0.70 and 0.89 at other monitoring points despite the fact that the UIB is a complex watershed to model, particularly on the daily time scale, as it covers a very vast area with a huge glaciated part and diverse hydro-climatic regimes as well as very sparse observation data. Table 4 indicates furthermore that, similar to the calibration period, the statistical results are also very good for the validation period. Table 4. Goodness-of-fit statistics for whole (calibration/validation) period.

Approach
The uncertainty in the UIB's future climate is evidently also reflected in the SWAT simulated projections of the future hydrology. However, as will be shown in the following subsections, a number of consistent patterns in the projected changes in UIB's future hydrology can be observed across the range of scenarios applied. To that avail, the calibrated SWAT model is forced with the downscaled future climate projections for the two RCPs (RCP 4.5 and RCP 8.5) of the five (5) selected CORDEX climate models (see Section 2.2.3). As listed in Table 2 earlier, each of the five (CORDEX-SA) RCMs selected represents essentially a particular future climate scenario, including the four possible future extreme scenarios groups (wet-warm, wet-cold, dry-warm dry-cold) for each RCP as well as of one scenario in each group, representing the mean future climate change over the UIB. For each scenario, the future hydrology is assessed for general flow characteristics, seasonality, high and low flows as well as extremes (return periods of extreme discharges) for two future periods, namely mid-century (2041-2070) and late-century (2071-2100). The results are summarized, along with other hydrometeorological variables and the glacier accumulation/melt (estimated with Equation (6)), both in absolute values and as relative changes compared with the historical reference period, in Table 5. The major features of this table are discussed in the subsequent subsections.

General Flow Characteristics
From Table 5, one can recognize that the projected changes in stream flow exhibit fairly consistent increases in the mid-as well as late-centuries for all future RCP/climate scenarios, except one, the dry-warm scenario (IPSL-MR-RCA4_RCP8.5 model), for which the simulation indicates a decrease in the mean daily flows (MDF) in both future periods. Apparently, in this dry-warm scenario, the increase in snow/glacier melt contributions to the flows (resulting from the obvious increase in temperature) is possibly balanced out and even surpassed by the decrease of snow accumulation (due to less snowfall owing to increased temperatures and, therefore, less snow melt) and increased losses from highly increased actual evapotranspiration (aET), which rises to 196 mm and 247 mm for the midand late-centuries, respectively, which for the latter represents more than a doubling from the aET value in the reference period. In fact, for this dry-warm scenario (IPSL_RCP 8.5_both periods), the possible increase in the simulated future flows owing to higher melting or precipitation seems to be overcompensated by as decrease owing to higher rates of aET (see shaded area in Table 5).

Response Surface Regression
The above findings regarding the relationships of both mean daily flow and aET with precipitation and temperature (Section 2.2.3) are explored further through "Response surface regression analysis/modelling (RSM)". Contour plots of the two response surfaces for these two variables are illustrated in Figures 5 and 6.
Regarding the "Flow" contour plot ( Figure 5), its pattern is in line with the earlier discussed quantitative flow results for the UIB in response to different climatic scenarios, namely a steady flow increase with increasing future precipitation, as long as the future increase in temperature is less than 6 • C. This is a clear indication of the dominance of the snow accumulation and snow/glacier melting process over the influence of losses due to increased aET. However, when the projected temperature changes are above 5 • C to 6 • C, the losses due to increased aET become more dominant, so that the flow decreases almost linearly with increasing temperature, and this holds for any amount of mean annual precipitation.
Similar behaviors are observable from the contour plot for "aET" (Figure 6) where, up to a projected temperature increase of 5 • to 6 • C, the change in aET is mainly affected by precipitation, whereas for temperature changes beyond the 5 • C mark, aET increases with a gradient of~30 mm/ • C so that, at 10-12 • C, the aET rate saturates at values around 350 mm/year. Interestingly, the changes of aET in the temperature range 5 • C-10 • C and the precipitation range 550-650 mm/year are similar for 1 • C changes in temperature and 100 mm/year changes in precipitation. In other words, in the named temperature range, an increase of aET by increasing the temperature by 1 • C is largely offset if the precipitation increases by 100 mm/year. Table 5. Future SWAT simulated hydrology (mean daily flows and low and high daily flow) and climate variables (aET, precipitation and temperature) and glacier melt in absolute numbers and relative to the reference period, obtained with climate drivers of the selected CORDEX GCM-Reg models under two RCPs for two future periods.

Approach
The uncertainty in the UIB's future climate is evidently also reflected in the SWAT simulated projections of the future hydrology. However, as will be shown in the following subsections, a number of consistent patterns in the projected changes in UIB's future hydrology can be observed across the range of scenarios applied. To that avail, the calibrated SWAT model is forced with the downscaled future climate projections for the two RCPs (RCP 4.5 and RCP 8.5) of the five (5) selected CORDEX climate models (see Section 2.2.3). As listed in Table 2 earlier, each of the five (CORDEX-SA) RCMs selected represents essentially a particular future climate scenario, including the four possible future extreme scenarios groups (wet-warm, wet-cold, dry-warm dry-cold) for each RCP as well as of one scenario in each group, representing the mean future climate change over the UIB. For each scenario, the future hydrology is assessed for general flow characteristics, seasonality, high and low flows as well as extremes (return periods of extreme discharges) for two future periods, namely mid-century (2041-2070) and late-century (2071-2100). The results are summarized, along with other hydrometeorological variables and the glacier accumulation/melt (estimated with Equation (6)), both in absolute values and as relative changes compared with the historical reference period, in Table 5. The major features of this table are discussed in the subsequent subsections.  In summary of this subsection, the generated contour maps of "Flow" (or "ET') can also be used as an easy first guess or prediction of flow in the UIB under diverse future precipitation and temperature scenarios, especially when advanced modeling and forecasting is not feasible. In summary of this subsection, the generated contour maps of "Flow" (or "ET') can also be used as an easy first guess or prediction of flow in the UIB under diverse future precipitation and temperature scenarios, especially when advanced modeling and forecasting is not feasible.

Monthly Flows and Seasonal Shifts
The UIB SWAT simulations show also changes in the future monthly and seasonal flows as indicated by the various panels of Figure 7. Thus, one may note changes in the yearly cycle and shifts in the timings of the high-and low-flow seasons, with similar trends but different magnitudes for the different scenarios.
For RCP 4.5 and all model/future period combinations considered, there is a shift of around one month in arrival of the peak flows, which is projected to be in June instead of July in the reference period. Simulations for RCP 8.5 show a similar trend, though with an even larger ahead-shift of up to two months. Also, for both RCPs and future periods, huge discharge increases are projected for the spring and early summer flows, so that the rise of the hydrograph in the future will start in spring rather than in midsummer, as for the reference period. The flow magnitudes in the spring season are also significantly higher than in the reference period (2 to 3 times for RCP 4.5 and 3-5 times for RCP 8.5). Furthermore, for all scenarios/periods, the future flow in autumn will be slightly increased. All these changes are most likely due to an increase in autumn or winter precipitation, leading to increased snow accumulation, and an earlier onset of snow/glacier-melt owing to rising temperatures.
For the RCP 8.5/late-century scenarios/periods, all five models project less increases and even a decrease for one model (IPSL) for the annual flows. The obvious reason is the higher rate of aET (Table 5), which negates the effect of increased precipitation and snow accumulation, so that for model IPSL which predicts less precipitation, the high aET induces considerable decreases in the average annual flows.
Our results above, for the intra-annual and magnitude changes of the flow, are in line with the projected hydrology of other studies (e.g., References [20,30,38]). Thus, a similar pattern of flows was reported by Lutz et al. [20] for RCP 8.5, namely, an earlier occurrence of higher-magnitude summer flows, particularly in the mid-century and less so towards the end of the 21st century.

Hydrological Extremes
Future changes in the hydrological extremes, i.e., low and high flows, are evaluated through their exceedance levels for different return periods (return period flow) as well as in terms of the 1st and 99th percentiles of the flows as listed previously in Table 5. For a better appreciation of the subsequent statistical frequency analysis, the hydrograph of the daily observed discharge at the UIB outlet station Bisham Qila is shown for the reference period  in Figure 8, with the various flow quantiles drawn in as horizontal lines. Noticeable are the extreme peaks. model IPSL which predicts less precipitation, the high aET induces considerable decreases in the average annual flows.
Our results above, for the intra-annual and magnitude changes of the flow, are in line with the projected hydrology of other studies (e.g., References [20,30,38]). Thus, a similar pattern of flows was reported by Lutz et al. [20] for RCP 8.5, namely, an earlier occurrence of higher-magnitude summer flows, particularly in the mid-century and less so towards the end of the 21st century.

Hydrological Extremes
Future changes in the hydrological extremes, i.e., low and high flows, are evaluated through their exceedance levels for different return periods (return period flow) as well as in terms of the 1st and 99th percentiles of the flows as listed previously in Table 5. For a better appreciation of the subsequent statistical frequency analysis, the hydrograph of the daily observed discharge at the UIB outlet station Bisham Qila is shown for the reference period  in Figure  8, with the various flow quantiles drawn in as horizontal lines. Noticeable are the The frequency histogram of the daily discharge is shown in Figure 9, which hints of a strongly skewed extreme value distribution. A similar histogram behaviour is found for the annual maximum series used for the subsequent analysis of the return periods. The frequency histogram of the daily discharge is shown in Figure 9, which hints of a strongly skewed extreme value distribution. A similar histogram behaviour is found for the annual maximum series used for the subsequent analysis of the return periods.

Annual Return Period Flows
For the evaluation of the annual return period (recurrence interval T) flow, exceedance probabilities (p = 1/T) of fitted Log-Pearson III value distributions (e.g., Reference [67]) to the observed (for the reference period) and the SWAT simulated future annual flow maxima using the various GCM Reg models' climate projections for the different RCPs/future periods were computed.
The four panels of Figure 10 show the annual return period flows obtained in this way, together with the ones for the reference period for the four RCP/future period combinations. Furthermore, the exceedance flows extracted for three return periods, i.e., 10, 15 and 30 years are listed in Table 6.
Both the figure and the table indicate remarkably more severe extreme discharges in the UIB compared to those in the reference period, in the mid-century and, depending on the model, often more so at the end of the century when exceedance flows are more than 100% higher than those in the reference period, particularly, for the more extreme RCP 8.5 scenario.
The different exceedance flow graphs in Figure 9 show similar patterns. Basically, for all tributaries and the main channel, the highest water discharge levels occur during the coinciding melting season (see Figure 7). Hence, the increases in the extreme return levels are essentially

Annual Return Period Flows
For the evaluation of the annual return period (recurrence interval T) flow, exceedance probabilities (p = 1/T) of fitted Log-Pearson III value distributions (e.g., Reference [67]) to the observed (for the reference period) and the SWAT simulated future annual flow maxima using the various GCM Reg models' climate projections for the different RCPs/future periods were computed.
The four panels of Figure 10 show the annual return period flows obtained in this way, together with the ones for the reference period for the four RCP/future period combinations. Furthermore, the exceedance flows extracted for three return periods, i.e., 10, 15 and 30 years are listed in Table 6.
Both the figure and the table indicate remarkably more severe extreme discharges in the UIB compared to those in the reference period, in the mid-century and, depending on the model, often more so at the end of the century when exceedance flows are more than 100% higher than those in the reference period, particularly, for the more extreme RCP 8.5 scenario.
The different exceedance flow graphs in Figure 9 show similar patterns. Basically, for all tributaries and the main channel, the highest water discharge levels occur during the coinciding melting season (see Figure 7). Hence, the increases in the extreme return levels are essentially regulated by the projected climate/hydrology changes during those months and are due to the higher snow/glacier melt owing to the projected temperatures increase and higher accumulation of snow during the preceding winter and spring while coinciding with intense rainfall events which can considerably exacerbate the peak flows (see column LDF 99 in Table 5).

Annual Return Period Flows
For the evaluation of the annual return period (recurrence interval T) flow, exceedance probabilities (p = 1/T) of fitted Log-Pearson III value distributions (e.g., Reference [67]) to the observed (for the reference period) and the SWAT simulated future annual flow maxima using the various GCM Reg models' climate projections for the different RCPs/future periods were computed.
The four panels of Figure 10 show the annual return period flows obtained in this way, together with the ones for the reference period for the four RCP/future period combinations. Furthermore, the exceedance flows extracted for three return periods, i.e., 10, 15 and 30 years are listed in Table 6.
Both the figure and the table indicate remarkably more severe extreme discharges in the UIB compared to those in the reference period, in the mid-century and, depending on the model, often more so at the end of the century when exceedance flows are more than 100% higher than those in the reference period, particularly, for the more extreme RCP 8.5 scenario.
The different exceedance flow graphs in Figure 9 show similar patterns. Basically, for all tributaries and the main channel, the highest water discharge levels occur during the coinciding melting season (see Figure 7). Hence, the increases in the extreme return levels are essentially regulated by the projected climate/hydrology changes during those months and are due to the higher snow/glacier melt owing to the projected temperatures increase and higher accumulation of snow during the preceding winter and spring while coinciding with intense rainfall events which can considerably exacerbate the peak flows (see column LDF99 in Table 5).  Table 6 indicates more quantitatively that the projected future peak flow magnitudes will increase between 50% and above 100% for the three annual return intervals investigated (10, 15 and 30 years). Similarly, for the even longer return periods of 50, 100 and 200 years, the peak flood amounts will more than double, and this holds for all climate scenarios/periods/models (see Figure  10). In other words, what has been a 30-, 50-, 100-or even 200-year flood event in the past reference period will become a flood with a return period of 10 years or below, while a 10-or 20-year flood will be witnessed every 2 years in the future.
This imminent future situation is very alarming because, even under the current flow regime, there are catastrophic flood events every 10 to 20 years, with catastrophic human and economic consequences. Thus, the projected increase in flood frequencies means that the UIB and, especially, its downstream basin are likely to witness more frequent and more intense flood events unless appropriate flood protection measures are put into place. Other studies (e.g., Reference [20]) have reported similar future trends in the projected extreme flows.   Table 6 indicates more quantitatively that the projected future peak flow magnitudes will increase between 50% and above 100% for the three annual return intervals investigated (10, 15 and 30 years). Similarly, for the even longer return periods of 50, 100 and 200 years, the peak flood amounts will more than double, and this holds for all climate scenarios/periods/models (see Figure 10). In other words, what has been a 30-, 50-, 100-or even 200-year flood event in the past reference period will become a flood with a return period of 10 years or below, while a 10-or 20-year flood will be witnessed every 2 years in the future.
This imminent future situation is very alarming because, even under the current flow regime, there are catastrophic flood events every 10 to 20 years, with catastrophic human and economic consequences. Thus, the projected increase in flood frequencies means that the UIB and, especially, its downstream basin are likely to witness more frequent and more intense flood events unless appropriate flood protection measures are put into place. Other studies (e.g., Reference [20]) have reported similar future trends in the projected extreme flows.  1 In the 30-year reference period. 2 Models full names and scenarios: CanESM2_RegCM4-4 (wet-warm), GFDL-ESM2M_RCA4 (wet-cold), IPSL-CM5A-MR_RCA4 (dry-warm), MPI-ESM-LR_RCA4 (dry-cold) and NorESM1-M_RCA4 (mean); 3 10-year recurrence interval, 4 15-year recurrence interval, 5 30-year recurrence interval.

High and Low Flows
For the analysis of various characteristics of the future low-and high-flow spells in the UIB, these spells have been defined here by the 1st (low daily flows (LDF 1 )) and 99th (high daily flows (DF 99 )) percentile historical flow discharges, as shown in Figure 8 and listed in Table 5, and which amount to threshold levels of 385 m 3 /s and 10,540 m 3 /s, respectively. In fact, the latter value corresponds roughly to a 2-year flood event for the reference period, as can be extracted from Figure 9.
The results of the SWAT simulations obtained with the various scenarios/periods/models are summarized in Table 7. From this table, one may notice similar future trends across all the assessed scenario combinations. More specifically, compared with the low-/high-flow values in the reference period (Table 7 and Figure 7), the magnitudes of mean low flow troughs at the 1st percentile have all decreased considerably (−13% to −19%), while the 99th percentile high flows all have significantly increased (10-25%) high spell peaks. Similarly, the magnitudes of the longest low flow spells for all the future scenarios have huge increases (21% to 132%), as do the magnitudes of the longest high flow spells, which have increased for all but the dry-warm scenario/model (RCP-8.5/IPSL-MR-RCA4), which predicts a huge future temperature increase but a precipitation decrease.
The results for the future high flow spells with discharges above the 99% threshold (10,540 m 3 /s) are in conformity with those of the previously computed projected return periods and flood magnitudes of floods (see Table 6), as the number of "high flow spells" have in most RCP/period/model cases increased by more than 300% in comparison to the reference period (Table 7). However, the mean duration of these future spells is considerably shorter by −15% to −69%, though some of them will also be much longer. All this means that these more frequent "high flow spell" events of shorter durations but with much higher peaks, i.e., flood concentration, will pose serious flood threats, as the present flood protection and retention infrastructures in the UIB are not designed for such extreme future floods.
In a similar manner, the 1% "low flow spells" (with troughs below 385 m 3 /s) are also projected to change drastically, as indicated by the appropriate columns of Table 7. Thus, the table shows that the UIB is expected to have more frequent and longer "low flow spell" events, especially for the wet-cold, dry-cold and mean model scenarios for both RCP/future periods. Although the wet-warm and dry-warm model scenarios will also have more frequent "low flow spell" events, these will have comparatively shorter durations than the others-possibly because of the high melt water contributions owing to the high future temperatures in these warm scenarios-but will still be much longer than those in the reference period.

Conclusions
Projecting future hydrology for the mountainous and highly glaciated UIB is a challenging task limited by many factors mainly due to uncertainties in the future climate projections and issues with the coverage and quality of available reference climatic data and hydrological modelling approaches. The current study has tried to investigate the impacts of future climate change on the hydrological behavior of the UIB during the 21st century while compensating for and avoiding, as much as possible, the limiting factors, so that the usability of the results could be enhanced as a viable set of future hydrological predictions.
The uncertainty in the future hydrological projections was addressed in two ways: (1) using reference data with improved quality, coverage and representation and (2) employing the wide range of different climate scenarios as provided by the various CORDEX GCM/RCM models in the SWAT model and getting diverse future hydrological projections. In addition, the full set of hydrological predictions was utilized to generate fitted Response Surface Regression Models (RSM) and contour maps for the two response variables "Mean Annual Flow" and "aET" as a function of the reference and projected mean annual precipitation and temperature. These RSMs can be used as a "hands-on" practical tool for preliminary forecast of the future hydrology in the UIB, under any temperature and precipitation projections, sourced from any climate model.
The overall results of the study are promising, since, notwithstanding the fact that there are always uncertainties in projecting future hydrological regime, almost all model/RCP scenarios have led to similar future trends of increased flows and intensified extremes. Thus, for all climate model scenarios under RCP 4.5, mean annual flows are projected to increase (relative to the 1976-2005 reference period) from 9.5-29.4% for the mid-century (2041-2070) period, while for the late century (2071-2100) period, a milder increase of only 10.2-19.1% is predicted. Similarly, all mean annual flow projections under RCP 8.5-except those based on model IPSL-MR-RCA4-show similar trends, with increases of 15.9-29.7% for the mid-and of only 9-27.5% for the late centuries. Here again, model IPSL-MR-RCA4 represents an exception, as its climate drivers lead to a steadily decreased SWAT simulated mean flow over the whole 21st century.
The seasonal shifts as well as the change in the extremes follow the same trends across all climate scenarios/models/period combinations analyzed. Hence, all projections show a future forward shift in the arrival of the peaks flows from, presently, July-August to late spring-early summer (May-June) as well as an increase of the magnitudes of the spring and winter flows.
The flow extremes are also projected to drastically intensify, with high (peak) floods having 50 to >100% more magnitude and tremendously decreased annual recurrence intervals, i.e., what has been something like a 100-year flood event or so in the past reference period will become a flood with a return period of 10 years or even below in the future. All this hint towards a tremendously increased flood hazard for the UIB in the future.
For the low flows, the more frequent and longer episodes with lower minimal flow values than those of the reference period are projected for all climate scenarios/models/period combinations, though these extreme alterations are somewhat milder for the future warm scenarios, most likely due to higher melt water contributions to the UIB flow discharge.