Prognosis of Hydro-Meteorological Attributes based on Simulation and Projection of Streamflow in a High-Altitude Basin using Hydrologiska Byråns Vattenbalansavdelning (HBV) Model

Subsistence of freshwater resources at high altitude regions has remained a paradox for stakeholder communities at both regional and global levels. To address such an issue, Hydrologiska Byrans Vattenbalansavdelning Light (hereafter HBV) model was used to assess hydro-meteorological shifts triggered under climate change scenarios in snow dominant region of Chitral river basin. The model performed well both during calibration and validation periods with Nash Sutcliffe Efficiency values of 0.91 and 0.81 respectively on daily time scale in the basin. The HBV was thereafter engaged for the projection of streamflow in the Chitral river basin using projected data of four statistically downscaled climate models with four emission scenarios for the 21st century. Multi-model ensemble projections of precipitation revealed an increase of up to 165% in monsoon inception period and an increase in temperature of up to 9.5°C in winter to summer transitioning period for the 2070‒2099 time slice under a high-end emission scenario. An increase of up to 122% in evapotranspiration was projected in the peak winter months for the 2070‒2099 time slice under the high-end emission scenario. Attributed to the significant increases in the temperature and the liquid precipitation, it was projected that basin streamflow had potential to increase by up to 182 % in the monsoon inception period for the 2070‒2099 time slice under the high-end emission scenario. It further indicated that precipitation might be falling as liquid rain most of the year, and snow will hardly accumulate in prognosticated future environements of the basin.


INTRODUCTION
The Hindukush-Himalayan (HKH) region hosts one of the prime assimilations of glaciers beyond Poles. They nourish several tributaries, which include Asia's paramount rivers -Kabul, Indus, Brahmaputra, and the Ganges. These waterways bluntly distress subsistence of 2.5 billion inhabitants of these basins. Himalayan snowmelt and glacial flux contribute at least 50% to the mean flow of significant watercourses of the region. For instance, during boreal winter, 70% of the flux from Kabul, Indus and Ganges rivers rests on the HKH meltwater. Towards the west side of China, glacial flux offers major freshwater supply during desiccated spell for up to 25% of inhabitants (Xu et al., 2007).
Environmental depreciation is one of the ultimate confrontations that the entire globe is encountering at present, particularly due to its influence on hydrological assets that can be disparate and indeterminate. In Pakistan where diverse microclimates exist, consequences of environmental degradation are essentially nontrivial. Pakistan's financial system centers along with crop production and thus is vastly reliant on the Indus irrigation system (Hewitt & Young, 1990). Fluctuations in stream flow levels are liable to upsurge conflicts among administrative units, specifically in downstream areas, concerning decreased flows in scorching periods and excessive flows followed by inundation in damp periods (Akhtar et al., 2008). Hence, in Pakistan, projected water estimation under greenhouse gas emission scenarios is indispensable for the design and control of hydrometeorological setups (Akhtar et al., 2008). Periodic streamflow prediction under various emission scenarios can deliver meaningful assistance for the supervision of national power schemes by presenting timely warning of excessive or deficient flow magnitudes which would necessitate offsetting of hydropower with thermal power sources (Archer & Fowler, 2008).
The hydrometeorology of rugged expanses in Pakistan is mainly reliant on the altitude of terrain, as highlands typically experience hefty sums of liquid and solid precipitation which may be deposited in form of glaciers. Massifs, likewise, signify exceptional sources for the discovery of climatic changes and estimation of climate-associated effects (Beniston, 2003). The intricacy and shared inter-reliance of high-altitude atmospheric routines on the subsistence of natural order create major complications in studying climate impacts (Beniston et al., 1997). Environmental evolution is projected to give feedback to augmented inconsistency in stream runoff owing to variations in state and strength of precipitation, in addition to the melting of deposited snow and permafrost. Runoff, at first, surges as snow thaws and then shrinks afterward when melting is halted (Fang & Pomeroy, 2007). A thorough exploration of snow and ice activities and their vulnerability to climate change has been done merely in partial extents of the HKH (Dahri et al., 2011). Limited findings in the HKH specify a sharp disparity in the comprehension of impacts borne due to climate change (Gautam et al., 2013). This apprehension gap exists in the HKH range of Pakistan too.
From the perspective of current changes in climate, knowledge of projected water resources is turning out to be beyond substance for policy devisers to earmark and assign volumes of water to provinces (Arnell, 2004;Shen et al., 2008). General Circulation Models (GCMs) have been perceived as the most intelligent implements for forecasting imminent consequences in the setting of climate evolution (Xu, 1999). Normally, the influence of climate evolution on the hydrometeorological scheme is assessed by expending GCMs as a feedback to an attuned hydrological model (Chiew, 2011;Kay et al., 2009;Teng et al, 2012;Wilby & Harris, 2006). In the present study, the HBV Bergstorm (1992) is applied to predict flow magnitudes in the Chitral river watershed under climate change settings recommended by Intergovernmental Panel on Climate Change (IPCC) established in Fifth Assessment Report (AR5). A limited number of hydrologists have deployed the HBV model for the estimation of projected watershed flow in both huge and minor basins. Akhtar et al. (2008) utilized HBV-Met (with downscaled climatic data) and HBV-PRECIS (with non-downscaled climatic data) to forecast potential streamflow magnitudes of Hunza river basin with diverse glacier extent scenarios. Sagar et al. (2017) acquired prognoses from the PRECIS-RCM and sourced them as feedback to the HBV model for estimation of imminent flow magnitudes of Karnali river watershed in Nepal for the period between 2030 and 2060 under A1B of Special Report on Emission Scenarios (SRES). Todorovic and Plavsic (2016) applied outputs of five GCM-RCM sequences, described by Langshot et al. (2013) as feedback in the HBV model to foresee impending flow magnitudes of Kolubara river basin, Serbia. Similarly, Usman et al. (2019) employed the HBV to simulate future streamflow in the Soan river basin of Pakistan.
The existing analysis connects the loopholes of preceding research works and offers a comprehensive investigation of the hydrological equilibrium of the Chitral river basin by deploying the HBV model (Konz & Seibert, 2010). Since the Chitral river basin constitutes about 13% of glaciers whose melted flux participates mostly in streamflow of the watershed, the preferred HBV model is practically fit for application over the glaciated basin. This paper supplements prior efforts (see e.g., Ali et al., 2015;Khalid et al., 2013;Lutz et al., 2016;Naeem et al., 2013;Shakir & Ehsan, 2016) and delivers expedient statistics to indigenous consultants for creating efficient strategies for impending sustainability and water regulations in the basin.
Primary intentions of this research are 1) to attune and authenticate the HBV model over the Chitral river watershed via in-situ hydro-meteorological records; 2) to project potential surge or decline of flows in the Chitral river by expending four statistically downscaled GCMs outputs under four Representative Concentration Pathways (RCPs) 2.6, 4.5, 6.0 and 8.5. This study has been carried out in the Chitral river basin, Pakistan between 1994 and 2099.

Description of the Study Area and Context
Chitral watershed is amongst the principal ancillary streams of the Indus watercourse that stems from immediate Baroghil pass in the Hindukush highlands from Kunhar tributary of Pakistan (Figure 1). The uppermost section of the Chitral watershed is known as Yarkhum river, and at downstream it is known as Chitral river located between 35°51ʹ48ʺ latitude and 71°47ʹ15ʺ longitude. Annual climatology of hydrometeorological parameters indicates highest precipitation of up to 2000 mm in winter retreating and spring approaching months, highest temperatures of up to 28°C in monsoon inception months, highest evapotranspiration of up to 5 mm in the monsoon inception months, and highest discharge of up to 882 m 3 /s over the basin (Figure 2). Efflux expanse of the Chitral watershed within Pakistan is 11,396 km². A normal spill of the watershed for 19 years 1994 -2012 is 108067 m 3 /s. Extreme of the measured water flux is 34601 m 3 /s recorded on the 6th of June 2005 and the minimum recorded discharge is 8 m 3 /s on the 12th of April 2004. The watershed is both ice and snow dispensed with stable fluxes round the annual cycle with supplementary overflow during warm and wet season i.e. July-September. It dispenses massive extents of snow resources in the Chitral basin and consequently throws in more than half of the discharge of the Kabul river basin.

Interventions made in the HBV
The current version of the HBV was retrieved from repositories of the Department of Geography, Zurich University (https://www.geo.uzh.ch/en/units/h2k/Services/ HBV-Model.html). Since the model could run at both daily and monthly time scales, we opted to simulate it on higher temporal resolution i.e. on a daily time scale. Input data required to drive the model was temperature and precipitation time series for every duration of increment, and climatology of both mean monthly temperature and potential evapotranspiration rates (Figure 3). It included an element that processed the feedback precipitation either as liquid rain or frozen snow, which relied on corresponding heat signature at every duration of increment. Any form of precipitation (in its presence) was subsequently treated in soil wetness unit wherever effective precipitation that added to the exterior spill  was calculated. The residual portion of the liquid precipitation contributed to the soil wetness accumulator which could subjectively be evaporated provided there was ample water substance below the surface. The HBV model's core output was water flux at the channel of the basin, with three sections, namely external flow, interflow (supply from immediate external flow) and baseflow (impact from subterranean flow) (Aghakouchak & Habib, 2010). The HBV was simulated under five main modules which are described as follows.

Data Collection
Hydrological data for the period 1994-2012 was obtained from the Surface Water Hydrology Project (SWHP) of Water and Power Development Authority (WAPDA, Pakistan), while meteorological data of Chitral observatory was provided by Pakistan Meteorological Department (PMD) for the same period. Evapotranspiration data was derived from (Irmak et al., 2003) using incoming solar radiation and mean temperature data provided by the PMD. An observation period of the meteorological data consisted of 19-years of historical daily weather records (1994-2012) ( Table 1). Data quality control was performed by identification of outlying (+/-3 standard deviations) values. The data were also checked to see if it were plausible over both time and space dimensions. Chitral station is located at 35°51ʹ48ʺ latitude and 71°47ʹ15ʺ longitude, with an elevation of 1497.8 meters above sea level (m.a.s.l.). Simulated and projected climate data of four statistically downscaled GCMs viz. GFDL-ESM-2M, MIROC-ESM-CHEM, IPSL CM5A-LR, and Nor-ESM1-M at 0.5 degrees horizontal and daily temporal resolution forced by RCPs 2.6, 4.5, 6.0, and 8.5 (Buda et al., 2016) was used to estimate future streamflow in the basin.

Snowmelt and Snow Accumulation Module
In this model routine, it was assumed that snowmelt and snow accumulation was directly proportional to temperature. Threshold temperature (TT) was taken as a significant model parameter -snow melted at temperatures above the TT and accumulated below it. Initially, the TT was set at 0°C (a reasonable assumption) and then was varied further within the parameter's sliding range. The accumulation of precipitation in the form of snow occurred if a precipitation event happened with temperature below the TT, otherwise, the input precipitation was assumed to be rainfall. The input precipitation did not contribute to a runoff if the temperature was below the TT. While in a scenario, where the temperature exceeded the threshold, both the snowmelt and the precipitation started contributing to the runoff. The snowmelt rate as water equivalent was estimated according to Eq. 1.
Where, is snowmelt rate as water equivalent / , is the degree-day factor /°/ , 1 is mean daily temperature °. As per definition, the indicates a reduction of rainwater substance in the snowpack triggered by a rise of 1°C beyond the diurnal ice limit.

Glacier Module
The degree-day method was used to compute glacial melt which was then added to the water content of the glacier. A small fraction Snow to Ice conversion factor (KSI) 1 of the simulated snow, storage was simulated to become glacial ice at each time step. Relationship of glacier melt reservoir out-flux varied diurnally to signify periodic evolution of subglacial outflow mechanism Eq. 2.
is the water equivalent of the snowpack on top of the glacier, and S is the water content of the glacier.

Effective Precipitation and Soil Vapor Sequencer
Rain or snow that fell upon the drainage basin was typically segregated into dual elements: the former contributed to permeation into the soil realm, and the later constituent added to external discharge. The module contributing to surface runoff was "effective precipitation" approximated by the HBV contingent on the soil dampness volume over a wet spell. Maximum soil moisture storage in the subsurface zone was described by field capacity ( ). The higher extent of soil vapor substance for liquid or solid precipitation enhanced the contribution of the precipitation to the runoff. As the soil vapor volume approached the , permeation decreased and the rainfall impact to the flow increased. The effective precipitation through an assimilation of the soil vapor substance was calculated using Eq. 3.
Where aliases effective precipitation , aliases net soil vapor , aliases limit of soil reservoir volume , aliases measure of each day's precipitation , and is a model parameter (shape coefficient) [-]. The controlled the amount of liquid water + which contributed to runoff. Moreover, as the approached the , the runoff coefficient increased exponentially.

Evapotranspiration Module
The long-term monthly mean potential evapotranspiration ( , = 1 − 12) was used as an input to calculate the actual evapotranspiration. Thereafter instead of respective time intervals in the replication cycle, corrected potential evapotranspiration was computed by reducing the hypothetical estimate centered on a deviation of daily mean temperature from normal of mean monthly temperature Eq. 4.
Where, is adjusted potential evapotranspiration , is a norm of daily temperature °, is normal of the mean monthly temperature of the °, and aliases model parameter

1°.
The soil vapor and the real evapotranspiration results were linked via soil Permanent Wilting Point (PWP). Eq. 5. expresses a balance amid soil vapor and real evapotranspiration.
Here is actual evapotranspiration . The Eq. 5 specifies that once the soil moisture equals the , the occurs at an equivalent pace as the EP . The aliases a soil vapor threshold intended for evapotranspiration i.e. whenever the is below the PWP, the remains below the EP . The volume of evapotranspiration reduces owing to a deficiency in soil vapor disposal lower than the according to the Eq. 5.

Runoff Response Module
Runoff was estimated at a watershed outlet based on a reservoir in this module. One reservoir was introduced to a model near-surface flow, whereas another was applied to simulate groundwater flux. From a sequential perspective, primary and second reservoirs simulated rapid and gentle subterranean mechanisms, respectively. The basins were exactly piped with one another via steady percolation pace . Once, fluid scale in the higher section of the basin exceeded its tolerance [ ], overspill followed swiftly from the higher section 0 . The flux reaction of additional dual channels was gentler. Control factors 0 , 1 , 2 responded to operations of the hierarchical reservoirs. Initial value of 0 was taken to be greater than 1 to ensure the fastest pace of the runoff. The response of the third outlet 2 was slower than the second outlet 1 and thus, 2 bore less value than 1 . The response functions of the outlets are given in Eq. 6.
is greater than aliases primary storage fluid rank , aliases secondary storage fluid rank , aliases limit of fluid rank .
The total simulated runoff , was attained by aggregating fluxes of the primary and the secondary sources; = ( 0 + 1 + 2 ).
Finally, the time-dependent model for a single linear reservoir was described as a catchment where runoff ( ) at time was supposed to be proportional to water storage ( ). The realization of a single linear reservoir was a basin with a porous outlet. Thus, we obtained a dependent and a scaled storage (or recession) coefficient equation from Darcy's Law Eq. 7.

Runoff Response Module
In the pre-rendering of the HBV, both the snowmelt and the ice flux and freeze routines were simulated. Additionally, to present a diversity of impacts from snow and ice fluxes, four-dimensional segregation in aspect genres comprising South, North, and East-West-Horizontal was additionally augmented ( Table 2). Glacier flux was emulated using analogous degree-day routine as of snow, nonetheless with the degree-day measure escalated for dissipation of ice in contrast to that of snow attributed to its depleted albedo.  prevalent glaciated expanse of about 1100 Km 2 in between 4500-5400 m.a.s.l. (Table 3). Optimum calibration values for the HBV parameters in the vegetation, the barren and the glaciation zones were retrieved by the Genetic Algorithm and Powell (GAP) optimization method (Seibert, 2000) and are presented in Table 4.
The coefficient of efficiency was determined using a statistical approach of the Nash-Sutcliffe model efficiency coefficient (NSE) which had been exercised widely to test prognostic control of hydrological emulators. It is defined as: Where is observed discharge, is simulated discharge and is the NSE. An value close to 1 indicates a virtually perfect fit. The HBV was found out to be very efficient both during calibration and verification with the NSE values of 0.91 and 0.81 for calibration and validation respectively over the basin (Figure 3). Table 3. Elevation wise area allocation of the three land cover zones  After obtaining optimum calibration and validation coefficients for the HBV, projected precipitation, temperature, and evapotranspiration based on the four GCMs and the four RCP scenarios were ingested in the HBV to project changes in the river flows for periods 2010-2039 (near), 2040-2069 (mid) and 2070-2099 (far) over the basin.

Projected Ensemble Changes in Precipitation
Projected changes in monthly precipitation of the engaged GCMs under the four RCPs for near, mid and far future time slices are shown in Figure 4. It is seen that the precipitation tends to increase in summers, and early autumn, with the IPSL-CM5A-LR suggesting an increase of up to 64 % in retreating monsoon periods for far-future time slice and a 49 % decrease in progressing monsoon periods for mid future time slice under the RCP 2.6 emission scenario. Under the RCP 4.5 emission scenario, the GFDL-ESM2M suggests an increase of up to 60 % in monsoon inception period precipitation in the near future time slice whilst the IPSL-CM5A-LR shows up to 54 % decrease in the autumn period in the far future time slice. An increase of 67 % in monsoon inception period for nearfuture time slice is projected according to the NoR-ESM1-M, while a decrease of 60 % in progressing monsoon period for far-future time slice is suggested by IPSL-CM5A-LR under the RCP 6.0 emission scenario. An extraordinary increase of 165 % in monsoon inception period for far-future time slice is projected by the GFDL-ESM2M, and a significant decrease of 87 % in peak winter periods for near-future time slice is projected according to MIROC-ESM-CHEM under the RCP 8.5 emission scenario. The overall consensus of models' ensemble project a significant increase in liquid precipitation (of the monsoon months) and stability in solid precipitation (of winter months) in the future time slices which is in-line with findings of Burhan et al. (2015) over the basin.

Projected Ensemble Changes in Temperature
Projected changes in monthly temperature under the four RCPs according to the engaged GCMs are presented in Figure  5. Temperatures have shown an overall increase till the end of the century in the Chitral river basin. According to the MIROC-ESM-CHEM, under the RCP 2.6 emission scenario, an increase of up to 4 °C is seen in winter to summer transitioning months for far-future time slice. As per the RCP 4.5 emission scenario, an even larger increase of 5.5 °C in the winter to summer transitioning months is projected according to the MIROC-ESM-CHEM results. Moreover, the RCP 6.0 emission scenario projects a 6.6 °C increase (which is greater than both the RCPs 2.6 and 4.5) in winter to summer transitioning months according to MIROC-ESM-CHEM for far-future time slice over the basin. Furthermore, a whopping 9.5 °C increase (highest increase amongst all the RCPs) is seen in the winter to summer transitioning period as suggested by MIROC-ESM-CHEM in the far future time slice projections under the RCP 8.5 emission scenario. The temperature projections of all the models are in a consensus of hypothesizing that the snow residence period for the metamorphic process of glaciation is susceptible to decrease owing to significantly high changes in all the emission scenarios over the basin. In other words, significant glacial melt is anticipated due to the inception of early summers projection in the basin (see e.g., Khalid et al, 2013).

Projected Ensemble Changes in Evapotranspiration
Projected changes in monthly evapotranspiration under the four RCPs according to the engaged GCMs are presented in Figure 6. A maximum increase of 59 % is seen in peak winter months of the mid and the far future time slices as displayed by the MIROC-ESM-CHEM, under the RCP 2.6 emission scenario. A higher increase of up to 79 % is seen in the peak winter months for the far future time slice as depicted by the MIROC-ESM-CHEM under the RCP 4.5 emission scenario. Additionally, even higher surge of up to 82 % is reflected in the peak winter months for the far future period as identified by the MIROC-ESM-CHEM, under the RCP 6.0 emission scenario. Furthermore, an outstanding increase of 122 % (which is highest amongst all the RCPs) is seen in the peak winter months for far-future time slice, as presented by the MIROC-ESM-CHEM, under the RCP 8.5 emission scenario. Projected increases in monthly evapotranspiration of the basin may be attributed to higher temperature changes in the Hindukush ranges as compared to those in the neighboring Himalayan ranges (see e.g., Lutz et al., 2016).

Projected Ensemble Changes in Stream Flow
In our analysis, the highest variability of flows is seen under the RCP 8.5 emission scenario. Projected changes in monthly discharge under the four RCPs according to the deployed GCMs are presented in Figure 7. An increase of up to 28 % in streamflow of monsoon progressing months is projected by the GFDL-ESM2M for the mid future time slice, while a decrease of up to 23 % in streamflow of peak winter months is projected by MIROC-ESM-CHEM for the far future time slice, under the RCP 2.6 emission scenario. Flows are seen to show an increase of up to 79 % in monsoon progressing months for a far-future time slice, and a decrease of up to 27 % in peak winter months for mid future time slice according to the IPSL-CM5A-LR under the RCP 4.5 emission scenario. Also, the IPSL-CM5A-LR projects an increase of up to 85 % in streamflow of monsoon progressing months, while it projects a decrease of up to 30 % in peak winter months for the far future time slice under the RCP 6.0 emission scenario. Moreover, an increase of up to 182 % is projected by the IPSL-CM5A-LR in monsoon inception period for a far-future time slice, and a decrease of up to 40 % is projected by the MIROC-ESM-CHEM for autumn flows under the RCP 8.5 emission scenario. It is to be noted that all GCMs based outputs agree on the direction of the flow changes under respective emission scenarios, however slight differences in magnitudes are seen in the stream flows for the projection periods. The results of the streamflow magnitudes depict an overall increase in summer and an overall decrease (stability) in winter flows. That may be attributed to a significant increase in temperatures of the spring leading to the monsoon months, as well as to the significantly high evapotranspiration changes in the winter months as per the GCMs based projections. However, the magnitude of increasing flows significantly overweighs the magnitude of the decreasing flows throughout the projections which is in line with findings of (Iqbal et al, 2018). Additionally, a significant increase in liquid precipitation (of the monsoon months) and stability in solid  (Bokhari et al, 2018) that too attributes to the direction of changes in the streamflow magnitudes for the projection periods.

Discussion on Possible Uncertainties
Although robust to a particular extent, the 0.5-degree resolution projections of possible climate changes may contain uncertainties for the moderately sized basin (only 11,396 km 2 ) located in the high altitude regions of Pakistan. This is illustrated in Figure 4 by variation in the estimates of changes in the normal of monthly total precipitation derived from different GCMs for each scenario and period. In current model settings, vaporization from the sea will transpire solely when there is no freeze. Glacier environments are designed with a straightforward scaling subroutine on atmospheric heat, which outturns in a delay connecting the atmospheric heat and sea temperature. Moreover, modeling of a catchment in high elevation zones might result in atypical snow accumulation and snow ablation in future projections. This might happen for areas above a certain elevation where, due to changes in projected lapse rate, temperatures are most of the year above the threshold temperature. Precipitation may be falling as liquid rain most of the year, and snow will hardly accumulate. This is a general issue in the currently engaged version of the HBV which does not take snowline displacement (e.g. due to changes in projected temperature) into account. The issues may lead to sources of uncertainty in the projections due to deficits of sea ice reduction and systematic increase of snow line in the HBV, which are foreseeable in case of the significant increase in the air temperature predicted by all the GCMs in their climatology.

CONCLUSION
Projected changes in the streamflow of the Chitral river basin are analyzed using the high resolution statistically downscaled GCMs output as input to the hydrological model HBV under the four AR5 based RCP emission scenarios. Chitral watershed is a sub-catchment of the greater Kabul river basin which is the second most vulnerable basin in South Asia after Figure 7. Projected climatology of discharge (mm/month) triggered with deployed GCMs under the AR5 based emission scenarios the mighty Indus river basin. Annual climatology of the hydrometeorological parameters indicates the highest precipitation of up to 2000 mm in the winter retreating and the spring approaching months, the highest temperatures of up to 28°C in the monsoon inception months, the highest evapotranspiration of up to 5 mm in the monsoon inception months, and the highest discharge of up to 882 m 3 /s over the basin. The HBV model engaged for projecting the streamflow comprises of the snowmelt and the snow accumulation module, the glacier module, the effective precipitation and the soil moisture module, the evapotranspiration module, and the runoff response module. The HBV is found out to be very efficient with the NSE values of 0.91 and 0.81 for calibration and validation respectively. Post-processed results project a significant increase in liquid precipitation (of the monsoon months) and stability in solid precipitation (of winter months) in the future time slices. The temperature projections of all the models are in a consensus of hypothesizing that the snow residence period for the metamorphic process of glaciation is susceptible to decrease owing to significantly high changes in all the emission scenarios over the basin. Moreover, the results also show that there is a significant projected increase in summer streamflow over the basin as per the models' output analysis. The results of this study are significant owing to emission based hydro-meteorological projections of the deficiently researched Chitral river basin and specifically to the multi-GCM based variation of outputs over the basin.

ACKNOWLEDGEMENT
The Authors are thankful to the PMD and the SWHP WAPDA for providing in-situ data for this study. Moreover, the authors are thankful to the SRTM for providing the DEM data through their repositories. Further, the authors are also grateful to USGS for sourcing remotely sensed indices data used in this study. Finally, the authors pay gratitude to Buda et al., (2016) for providing statistically downscaled projections of the GCMs used in this study.