Simulating Current and Future River-Flows in the Karakoram and Himalayan Regions of Pakistan Using Snowmelt-Runo ff Model and RCP Scenarios

Upper Indus Basin (UIB) supplies more than 70% flow to the downstream agricultural areas during summer due to the melting of snow and glacial ice. The estimation of the stream flow under future climatic projections is a pre-requisite to manage water resources properly. This study focused on the simulation of snowmelt-runoff using Snowmelt-Runoff Model (SRM) under the current and future Representative Concentration Pathways (RCP 2.6, 4.5 and 8.5) climate scenarios in the two main tributaries of the UIB namely the Astore and the Hunza River basins. Remote sensing data from Advanced Land Observation Satellite (ALOS) and Moderate Resolution Imaging Spectroradiometer (MODIS) along with in-situ hydro-climatic data was used as input to the SRM. Basin-wide and zone-wise approaches were used in the SRM. For the zone-wise approach, basin areas were sliced into five elevation zones and the mean temperature for the zones with no weather stations was estimated using a lapse rate value of −0.48 ◦C to −0.76 ◦C/100 m in both studied basins. Zonal snow cover was estimated for each zone by reclassifying the MODIS snow maps according to the zonal boundaries. SRM was calibrated over 2000–2001 and validated over the 2002–2004 data period. The results implied that the SRM simulated the river flow efficiently with Nash-Sutcliffe model efficiency coefficient of 0.90 (0.86) and 0.86 (0.86) for the basin-wide (zone-wise) approach in the Astore and Hunza River Basins, respectively, over the entire simulation period. Mean annual discharge was projected to increase by 11–58% and 14–90% in the Astore and Hunza River Basins, respectively, under all the RCP midand late-21st-century scenarios. Mean summer discharge was projected to increase between 10–60% under all the RCP scenarios of midand late-21st century in the Astore and Hunza basins. This study suggests that the water resources of Pakistan should be managed properly to lessen the damage to human lives, agriculture, and economy posed by expected future floods as indicated by the climatic projections.


Introduction
The change in climate patterns is one of the major challenges that the world is confronting today, especially due to its diverse and uncertain impacts on water resources.Anthropogenic warming has caused the retreat of glaciers worldwide since the industrial revolution.An effect of global warming is significantly attributed to the increase in the glacial retreat rate [1][2][3] across the world except for the Karakoram glaciers, which are reported to be advancing [4][5][6][7] in a few regions such as northern areas of Pakistan.The impact on access of water resources for irrigation, drinking, hydro-based industries, and hydropower are the major implications arising from changes in peak runoff and water balance [8,9] due to climate change.Climate variability is likely to influence the accessibility and intensity for irrigation and hydroelectric power in a densely-populated country like Pakistan.There is already a conflict between the provinces of Pakistan, resulting from inaccessibility and inefficient supply of water, and the potential fluctuations in water availability can escalate the dispute.Indus is the major river to provide with most of the water supplies in downstream areas of Pakistan.Catchment area of Indus upstream Tarbela dam (Figure 1) is called the Upper Indus Basin (UIB) residing in the Hindukush-Karakoram-Himalayan (HKH) region and abstracts most of the water supply as a result of snow and glacier melting.An authentic assessment of future water availability under projected climate change from UIB is a pre-requisite for devising, management, and operation of various hydrological infrastructures in Pakistan [10][11][12].
Water 2019, 11, x FOR PEER REVIEW 4 of 19 flows.Due to the increased runoff from UIB, Pakistan observed consecutive flooding from 2010 to 2016, affecting millions of people, infrastructure, and economy.UIB is characterized by large altitudinal variations over small horizontal distances and has a mean elevation of about 4750 m ASL [45].The region, particularly above 4000 m ASL generating most of the flow to UIB, has a dearth of in-situ climate data [54,61], primarily due to lack of observatories.Even if there are few observatories, they are installed at valley floors and the highest climate station is located at 4440 m ASL in Khunjerab (Figure 1) [45].Secondly, accessibility of data is difficult due to the rugged terrain and remote location of the region.The study sites consist of two comparable sub-basins in the UIB; the Astore basin (Western Himalaya) and the Hunza basin (Karakoram Range) as shown in Figure 1.SRM ability to simulate the current and future snowmelt-runoff from these two basins was evaluated in this study.These two sub-basins were selected because of their different geographical location within the bigger UIB (Figure 1) and slightly different climates such as Astore receives more amount of precipitation than the Hunza basin as observed from the available data records.Both basins do not have any weather station above mean elevations so more than 50% of their area remains un-gauged so the true amount of precipitation, especially at higher altitudes is unknown.Location of studied catchments, Digital Elevation Model (DEM) showing the elevation zones of the study area and available gauging stations (flow gauge and weather stations) are presented in Figure 1.Astore basin covers an area of 3990 km 2  as estimated from the DEM.The mean elevation is about 4100 m ASL as determined using the hypsometric curve (Figure 2) with maximum altitude reaching above 8000 m ASL (Figure 1).Approximately, 6% (245 km 2 ) of the catchment's area is glacier covered [26].Hunza basin has an area Snow cover is one of the most significant variables that can be accustomed to water resource management [13,14].Snow cover areas need to be measured properly in order to estimate and forecast the runoff from regions [15] such as UIB, but there is a paucity of comprehensive and in-detail study of ice and snow processes and their relationship with the hydrological regime in this region.In-situ data is limited to lower and accessible altitudes and may not truly represent the snow and glacier changes [4] at higher altitudes like those of UIB.Remote sensing data and various modeling approaches such as physical-based models [13], cellular automata (CA) models [16,17], and artificial neural networks [18,19] were used previously to analyze the snow cover in different regions of the world.Remote sensing data offers the quantitative examination of physical properties of snow and glaciers in remote areas where accessibility of data is expensive and dangerous [20].Estimates based on remote sensing are indirect and require validation [21].Validation of remotely sensed data is often difficult in the absence of in-situ data; therefore, the use of multisource observation is a common practice for evaluating the validity of remotely sensed data [22].Satellite snow cover maps of the National Oceanic and Atmospheric Administration (NOAA) provide weekly data that shows absence or presence of snow on a grid [14].Many studies used Moderate Resolution Imaging Spectroradiometer (MODIS) daily, 8-day, and monthly snow cover products produced by the National Snow and Ice Data Center (NSIDC) for snow cover estimation [22][23][24][25][26][27][28].It is considered a real standard for snow cover estimation.Impacts of future climate projections on the snow cover area were studied using the cellular automata (CA) model [16], supported with MODIS satellite data.Daily snow products produced by MODIS sensors were available since 2000.More than 90% efficiency of MODIS products is documented by scientists [29][30][31].These studies endorse the application of MODIS data for determining snow cover area at high elevations, and in the mountainous area of UIB.There is a suite of MODIS snow products which are made available for free by the NSIDC and NASA's online portal.
The remotely sensed snow cover data (such as MODIS) may be used as a key input to simulate the current and future snowmelt runoff from mountainous regions like UIB.Change in glaciers-and snowmelt-runoff under future climate change may adversely affect the water supplies of downstream dependent populations [15].It is, therefore, necessary to model the influence of climate variability on the snow and glacier melt runoff.Many models (e.g., SSARR, HEC-1, NWSRFS, PRMS, SRM and GAWSER etc.) including Snowmelt-Runoff Model (SRM) were used in the mountainous regions of the world [32] to simulate the snowmelt runoff for operational uses.SRM was developed by Martinec [33] to simulate the stream flow in mountainous basins where meltwater from snow is the primary runoff source.SRM uses meteorological data, and the remotely sensed snow cover map as input variables.The simulation efficiency of the SRM model was less responsive to the daily precipitation data because it was designed to be used for remote mountainous areas having scarcity and well-known gauging errors for in-situ precipitation measurements [34,35] at high altitudes.As compensation to the lack of precipitation data, SRM computes discharges using snow-cover data as a key input [36].It was used extensively in the poorly gauged catchments [27] due to its low data requirements, simplicity, and good performance.This model was used successfully in more than one hundred mountainous basins [35], including those of HKH region [28,[36][37][38][39]. SRM was tested by World Metrological Organization [32] and it was used to investigate the impacts of climatic changes on seasonal river flows in North America, the Swiss Alps, the Himalaya, and in the Upper Indus Basin [40,41].
Climate varies greatly within the UIB due to a large catchment area and altitude difference.UIB is influenced by two precipitation regimes i.e., winter precipitation brought by westerly circulation in Hindukush and Karakoram regions, and the monsoon rainfall (summers) in the Himalayan region of the basin [42,43].The eastern UIB receives more monsoon precipitation as compared to the western UIB; therefore, the northern and western parts have distinct climate from that of the eastern parts of UIB [44,45].Therefore, estimating snowmelt runoff at sub-basin and altitudinal zones scale may be more meaningful as compared to the entire basin.SRM offers the flexibility to study snowmelt runoff at small basin (basin-wide) and elevation band (zone-wise) scales.
To study the impacts of climate change on the snowmelt-runoff, the climatic projections generated by global and regional climate models may be coupled with SRM.For example, the fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC) introduced a new set of greenhouse gas concentration climate change scenario projections [46] recognized as Representative Concentration Pathways (RCPs).The scenarios were projected for a near term up to 2035 and long term up to 2100.There were four scenarios which were developed and named after the level of radiative forcing, measured in Watts per square meter (W/m 2 ) by the year 2100.These scenarios contained a mitigation scenario (RCP 2.6), two medium stabilization scenarios (RCP 4.5 and RCP 6), and a high emission scenario (RCP 8.5).The RCP 8.5 appeared from little effort to lessen emissions and indicate a failure to limit warming by 2100 [47].The output of physical models could be treated using statistical techniques to generate future climate scenarios [48].During the simulation of current climate, the global and regional climate models showed many types of biases [49].Subsequently, the model-projected future climate seldom provided an authentic simulation of the future climate.To lessen the biases of models on the simulations, different methods were used while developing the climate projections such as bias correction and delta change approaches [50].
The assessment of potential impacts of climate change on snowmelt-runoff using the statistically downscaled scenarios for Indus basin are not yet studied in the sub-basins of UIB.Therefore, present study focused on: (1) evaluating the ability of snowmelt runoff model to simulate river flow in the sub-basins of UIB under different climate regimes, and (2) studying the impacts of future climatic changes on the stream flow of these sub-basins using statistically downscaled RCP scenarios, for an efficient water management and planning system in the country.

Case Study
The economy of Pakistan is dependent on the agriculture sector, which thrives on the water furnished by the mighty Indus River and its tributaries [51].The Indus basin having an area of ≈1.1 million km 2 originates from the Tibetan Plateau and the HKH region and extends across Pakistan (47%), India (39%), China (8%), and Afghanistan (6%) with its upper portion (called UIB) resting upstream of Tarbela reservoir [52].The total catchment area of the UIB is ≈200,000 km 2 and almost its 12% area is covered with glacial ice [45].More than 60% flow concentrated from Indus River at Tarbela is supplied from the melting snow and glacial ice [8,45,[53][54][55][56].The sub-basins of the UIB are the Gilgit, Astore, Hunza, Shingo, Shyok, Shigar, and Zanskar basins.UIB is sustained by glacier melt, seasonal snowmelt, and direct runoff from winter and summer monsoon rainfall [38,52,57].The lower Indus Basin is arid and semi-arid and is dependent on the water from UIB [53].Approximately 800 million people living downstream rely on supplies by the UIB for their freshwater needs [58,59].
Cryosphere melt contributions are significant (>60%) in the UIB, and therefore it is highly susceptible to climate change.An increase of 1 • C in the temperature would raise the summer runoff by 16-17% in the UIB [60] as a result of snow and glacial melt.Continued melting of the glacial ice may result in severe flooding in this region and shrinking glaciers will be followed by fewer river flows.Due to the increased runoff from UIB, Pakistan observed consecutive flooding from 2010 to 2016, affecting millions of people, infrastructure, and economy.UIB is characterized by large altitudinal variations over small horizontal distances and has a mean elevation of about 4750 m ASL [45].The region, particularly above 4000 m ASL generating most of the flow to UIB, has a dearth of in-situ climate data [54,61], primarily due to lack of observatories.Even if there are few observatories, they are installed at valley floors and the highest climate station is located at 4440 m ASL in Khunjerab (Figure 1) [45].Secondly, accessibility of data is difficult due to the rugged terrain and remote location of the region.
The study sites consist of two comparable sub-basins in the UIB; the Astore basin (Western Himalaya) and the Hunza basin (Karakoram Range) as shown in Figure 1.SRM ability to simulate the current and future snowmelt-runoff from these two basins was evaluated in this study.These two sub-basins were selected because of their different geographical location within the bigger UIB (Figure 1) and slightly different climates such as Astore receives more amount of precipitation than the Hunza basin as observed from the available data records.Both basins do not have any weather station above mean elevations so more than 50% of their area remains un-gauged so the true amount of precipitation, especially at higher altitudes is unknown.Location of studied catchments, Digital Elevation Model (DEM) showing the elevation zones of the study area and available gauging stations (flow gauge and weather stations) are presented in Figure 1.Astore basin covers an area of 3990 km 2 as estimated from the DEM.The mean elevation is about 4100 m ASL as determined using the hypsometric curve (Figure 2) with maximum altitude reaching above 8000 m ASL (Figure 1).Approximately, 6% (245 km 2 ) of the catchment's area is glacier covered [26].Hunza basin has an area of about 13,733 km 2 and nearly 25% (3417 km 2 ) area is covered with glaciers [26].Hunza basin has a mean elevation of ≈4650 m ASL which was determined using the hypsometric curve (Figure 2) estimated from the DEM.Mean monthly temperatures from 2000-2013 of Astore and Hunza basins are presented in Figure 3a.Mean monthly minimum temperature for the Astore (−1.9 °C) and Hunza (−11.7 °C) basins was observed in January while mean monthly maximum temperature was observed in July (11.5 °C) and August (20.7 °C) for Hunza and Astore, respectively (Figure 3a). Figure 3b shows the monthly precipitation totals recorded in the Astore and Hunza River basins from 2000 to 2013.The Figure 3b shows that maximum precipitation was received during the September-April, in winter and spring season, as snow, because the temperature during these months was below or near freezing (Figure 3a) in both basins.Most of this precipitation originates from Westerlies.
Astore and Hunza rivers are gauged at Doyian and Dainyor hydrometric stations (Figure 1), respectively.Astore river joins the Indus just downstream of its flow gauge.The Hunza river joins the Indus River few kilometers downstream after the confluence with the Gilgit River at Dainyor.River discharges converted to equivalent water depths (in mm) upon dividing by the basin area for both basins are presented in Figure 3c, which represents that Astore has more runoff than the Hunza which is also indicative of the more precipitation amounts (Figure 3b) in the Astore than Hunza basin.Astore also has higher temperatures than the Hunza (Figure 3a), so it may have more melt runoff generation than the Hunza basin.Astore River flow starts increasing in the month of March, and a peak discharge is achieved in June, as observed over the available data period (Figure 3c) when monsoon rainfall is superimposed by the peak melting season.The Hunza River flow starts increasing in the month of April and reaches its maximum during July as observed over the available data period (Figure 3c).The reason for delayed flow rising limb in the Hunza River (in April as compared to March in Astore river) is the mean temperature which only starts to reach above freezing in the month of April (Figure 3a).Mean snow cover area estimated using the MODIS snow product for the period of 2000-2016 is presented in Figure 3d.Hunza basin has slightly higher mean annual SCA (64%) than the Astore basin (57%).Maximum mean snow cover area is found to be in the March (end of snow accumulation season) and minimum in August (end of snowmelt season) in these basins (Figure 3d).3a). Figure 3b shows the monthly precipitation totals recorded in the Astore and Hunza River basins from 2000 to 2013.The Figure 3b shows that maximum precipitation was received during the September-April, in winter and spring season, as snow, because the temperature during these months was below or near freezing (Figure 3a) in both basins.Most of this precipitation originates from Westerlies.
Astore and Hunza rivers are gauged at Doyian and Dainyor hydrometric stations (Figure 1), respectively.Astore river joins the Indus just downstream of its flow gauge.The Hunza river joins the Indus River few kilometers downstream after the confluence with the Gilgit River at Dainyor.River discharges converted to equivalent water depths (in mm) upon dividing by the basin area for both basins are presented in Figure 3c, which represents that Astore has more runoff than the Hunza which is also indicative of the more precipitation amounts (Figure 3b) in the Astore than Hunza basin.Astore also has higher temperatures than the Hunza (Figure 3a), so it may have more melt runoff generation than the Hunza basin.Astore River flow starts increasing in the month of March, and a peak discharge is achieved in June, as observed over the available data period (Figure 3c) when monsoon rainfall is superimposed by the peak melting season.The Hunza River flow starts increasing in the month of April and reaches its maximum during July as observed over the available data period (Figure 3c).The reason for delayed flow rising limb in the Hunza River (in April as compared to March in Astore river) is the mean temperature which only starts to reach above freezing in the month of April (Figure 3a).Mean snow cover area estimated using the MODIS snow product for the period of 2000-2016 is presented in Figure 3d.Hunza basin has slightly higher mean annual SCA (64%) than the Astore basin (57%).Maximum mean snow cover area is found to be in the March (end of snow accumulation season) and minimum in August (end of snowmelt season) in these basins (Figure 3d).

Datasets
Two types of data sets were used to simulate the stream flow using SRM in the Astore and Hunza basins.These datasets consist of remotely sensed satellite data, including DEM and snow cover maps, and the field observation records consisting of hydro-meteorological variables such as mean temperature, precipitation, and the stream flows.
Two types of remote sensing data were used: (a) DEM data of Advanced Land Observation Satellite (ALOS); and (b) MODIS (MOD10A2) snow product.ALOS DEM was used to delineate the catchment and zonal areas, and to study the topographical characteristics such as elevation variation and hypsometric curves of the study area.ALOS has 12.50-m spatial resolution and was downloaded from Alaska Satellite Facility's online data portal (https://vertex.daac.asf.alaska.edu/).MODIS 8-daily snow product (MOD10A2) was acquired from the Terra satellite operated by NASA.This data was downloaded from the online data portal of the National Aeronautics and Space Administration (NASA) and the United States Geological Survey (USGS).This snow cover product was validated, and its suitable accuracy was reported by many studies [23][24][25]28,30].MOD10A2 have fields for maximum snow cover area over 8-day duration at a spatial resolution of nearly 500 m [29].Data available from the year 2000 to 2016 was downloaded for the study area.
Daily climatic data was obtained from the Water and Power Development Authority (WAPDA), and Pakistan Meteorological Department (PMD) of Pakistan.There are three weather stations in Astore (Rama, Rattu, and Astore) and three in Hunza (Khunjerab, Naltar, and Ziarat).Long timeseries (≈50 years) is available for PMD stations and a medium time-series (≈20 years) is available from WAPDA stations.Measurements of stream discharge, in Pakistan are recorded by the SWHP (Surface Water Hydrology Project) of WAPDA with the earliest observations dating from 1960 [45,62].Daily

Datasets
Two types of data sets were used to simulate the stream flow using SRM in the Astore and Hunza basins.These datasets consist of remotely sensed satellite data, including DEM and snow cover maps, and the field observation records consisting of hydro-meteorological variables such as mean temperature, precipitation, and the stream flows.
Two types of remote sensing data were used: (a) DEM data of Advanced Land Observation Satellite (ALOS); and (b) MODIS (MOD10A2) snow product.ALOS DEM was used to delineate the catchment and zonal areas, and to study the topographical characteristics such as elevation variation and hypsometric curves of the study area.ALOS has 12.50-m spatial resolution and was downloaded from Alaska Satellite Facility's online data portal (https://vertex.daac.asf.alaska.edu/).MODIS 8-daily snow product (MOD10A2) was acquired from the Terra satellite operated by NASA.This data was downloaded from the online data portal of the National Aeronautics and Space Administration (NASA) and the United States Geological Survey (USGS).This snow cover product was validated, and its suitable accuracy was reported by many studies [23][24][25]28,30].MOD10A2 have fields for maximum snow cover area over 8-day duration at a spatial resolution of nearly 500 m [29].Data available from the year 2000 to 2016 was downloaded for the study area.
Daily climatic data was obtained from the Water and Power Development Authority (WAPDA), and Pakistan Meteorological Department (PMD) of Pakistan.There are three weather stations in Astore (Rama, Rattu, and Astore) and three in Hunza (Khunjerab, Naltar, and Ziarat).Long time-series (≈50 years) is available for PMD stations and a medium time-series (≈20 years) is available from WAPDA stations.Measurements of stream discharge, in Pakistan are recorded by the SWHP (Surface Water Hydrology Project) of WAPDA with the earliest observations dating from 1960 [45,62].Daily stream flow data of the Astore and Hunza rivers were obtained from the WAPDA for the last 30 years.

Methods
SRM was applied to simulate the current and forecast the future river flows from the Astore and Hunza basins under climate variability scenarios in this study.The parameters and variables used in SRM are presented in Table 1.SRM uses daily snow cover, temperature, precipitation, and catchment characteristics as input variables.The model was run over the five-year data period of 2000-2004.This period was chosen because we had a continuous available data for this time period for both the catchments for each of the variables mentioned in Table 1 including daily river discharge.The basic parameters of the model were calibrated over 2000-2001 period and validated over 2002-2004 to simulate the river flows efficiently on both studied river basins.Model efficiency was estimated using Nash-Sutcliff (NS) model efficiency coefficient and the difference of volume (D v ) in the measured and simulated river flow.The calibrated and validated model was then used to investigate the influence of climate change on the river flow from the study area under the standard climate variability scenarios published by IPCC and statistically downscaled by [49] for the Indus basin.
The RCP scenarios for the mid-and late-21st century, downscaled statistically [49] for the Indus River Basin used in this study, were corrected for biases.Simulations of Coupled Model Inter-comparison Project phase 5 (CMIP5) for the Indus Basin were compared with the datasets of Climate Research Unit (CRU) and Asian Precipitation-Highly-Resolved Observational Data Integration Towards Evaluation (APHRODITE) [49].While downscaling the RCP scenarios for the Indus Basin, the biases were corrected systematically between simulated and observed datasets using the Equidistant Cumulative Distribution Functions Matching method (EDCDFm) and the high-resolution simulations were then statistically downscaled [49].Afterwards, the temperature and precipitation were projected for the mid-and late-21st century for the Indus basin [49].The change in mean annual temperature and precipitation under the RCP scenarios for the Indus basin used in this study is shown in Table 2.It is recommended to apply SRM by dividing the catchment area into zones (zone-wise) of 500 m elevations if the altitude varies largely in the catchment [35].SRM was applied for entire basin (basin-wide) and for zones (zone-wise) in both catchment areas in this study.Most of the studies used both methods for simulation considering the fact that the physical environment changes with altitude [27,34,36].
For the basin-wide simulation of the river flow, the mean temperature and total precipitation were calculated from the available daily data of climate stations.The daily snow cover was calculated using linear interpolation from the 8-day snow cover obtained from MODIS MOD10A2 images.Linear interpolation was previously used by many studies for extracting daily snow cover area from the MODIS dataset [13,27,36,63].In addition to climatic data, basin area and mean elevation estimated from the DEM were also used as inputs to the SRM.For the zone-wise application, both basins (Astore and Hunza) were divided into five zones according to elevation, with a difference of about 1000 m in the mean elevation (hypsometric) between every zone.Physical characteristics of these elevation zones used in the SRM as input for both basins are presented in Table 3.Average of daily precipitation totals estimated at weather stations of each basin (as estimated for the basin-wide application) was used in each zone because ground data was not available from each zone.SRM is less sensitive to precipitation which is often under-estimated or scarcely gauged in the high-altitude regions.It rather takes SCA and mean temperature as key inputs to simulate the snow melt generated runoff.Daily mean temperature data available from the weather stations situated inside the basins was extrapolated using an estimated lapse rate value to get the mean temperature values for the zones where there is no weather station.Lapse rate value was estimated using the available data and elevation of installed weather stations in the study area.A lapse rate value range of −0.48 • C to −0.76 • C/100 m was found between different elevation zones and weather stations.To obtain the zonal snow cover for the zone-wise application of the SRM, MODIS snow images were reclassified according to the zones presented in Table 3.The snow covered area was then estimated from the 8-daily snow images and linearly interpolated to get daily SCA for each zone throughout the time period of data, for the Astore and Hunza basins.This zonal SCA is a major input to simulate the snowmelt runoff using zone-wise approach of SRM.Zonal mean temperature extrapolated from available ground stations using a temperature lapse rate value and snow-covered areas estimated from reclassified MODIS snow images for zone-wise approach are presented in the Figures S1 and S2.

Results and Discussion
The results acquired after applying SRM to simulate the Astore and Hunza River flow for basin-wide and zone-wise approaches are presented in the Sections 4.1 and 4.2.The model was calibrated on two years (2000)(2001) and validated using the data of three years (2002)(2003)(2004) as explained in the Section 3.2.The calibrated model was then employed to study the future influence of climate change on the river discharge from Astore and Hunza basins and results are presented in Section 4.3.

Simulation of the River Discharge by SRM Using the Basin-Wide Approach
SRM was calibrated over the 2000-2001 period to estimate the best suited parameters to simulate the basin-wide river flow from Astore and Hunza basins.The values of these parameters were applied to validate the model using the data of three years (i.e., [2002][2003][2004].The estimated parameters for Water 2019, 11, 761 9 of 19 basin-wide application of SRM for the study area are presented in Table 4.The values of parameters were either estimated while calibration of SRM or taken from the previously published literature for this region.For example, range of parametric values for Degree-Day Factor, Critical temperature and Lag time were extracted from previous literature and rest were either estimated from the available data or during the calibration process.The efficiency of the SRM to simulate daily flow with basin-wide application in the Astore and Hunza River basins for the 5 years (2000)(2001)(2002)(2003)(2004) in terms of statistical analysis (Nash-Sutcliffe coefficient, difference of volume, Pearson and Kendall's rank correlation tests) is presented in Table 5.For the Astore River basin, the most efficient application of the SRM model was for the calibration year 2000 with the NS coefficient of 0.95 and −0.06% difference in volume (Table 5).NS value was 0.82 or higher which shows a satisfactory efficiency of SRM in the Astore basin.The correlation coefficient values were not less than 0.82 and 0.68, for the Pearson and Kendall rank correlation tests, respectively, in the Astore basin.
In the case of Hunza River basin, the SRM model showed best efficiency for the calibration year 2000, with 0.97 NS coefficient value and −0.76% difference in volume (Table 5).The value of NS was not below 0.73 which shows a good efficiency of SRM, but this value is less than that found for the Astore basin (0.82).The Pearson and Kendall's rank correlation tests were not less than 0.87 and 0.80, respectively (Table 5).
Overall graphic representation of the basin-wide SRM efficiency during calibration and validation period for the Astore basin is presented in Figure 4a.The NS coefficient for the calibration and validation period was 0.93 and 0.87, respectively (Figure 4a).Measured and simulated flow curves indicate that SRM simulated the high and low flows very efficiently during almost the entire data period.Value of Pearson coefficient was 0.95 for the correlation between the measured and simulated flow over the entire data period (Figure 5a) in the Astore basin.Scatter plot (Figure 5a) indicated a very good agreement between the low, mid, and high level measured and simulated flow in the Astore basin.Overall graphic representation of the basin-wide SRM efficiency during calibration and validation period for the Astore basin is presented in Figure 4a.The NS coefficient for the calibration and validation period was 0.93 and 0.87, respectively (Figure 4a).Measured and simulated flow curves indicate that SRM simulated the high and low flows very efficiently during almost the entire data period.Value of Pearson coefficient was 0.95 for the correlation between the measured and simulated flow over the entire data period (Figure 5a) in the Astore basin.Scatter plot (Figure 5a) indicated a very good agreement between the low, mid, and high level measured and simulated flow in the Astore basin.
Graphic representation of the basin-wide SRM efficiency during calibration and validation period for the Hunza basin is presented in Figure 4b.The value of NS coefficient for the calibration and validation period was found to be 0.86 (Figure 4b).The measured and simulated flow curves indicated that SRM simulation at high and low flows was efficient during the whole study period.Pearson coefficient value was 0.92 for the correlation between the measured and simulated flow over the entire data period (Figure 5b) in the Hunza basin.It was also noticed that SRM overestimated the mid-level flows and underestimated the high flows in the Hunza basin (Figure 5b).Graphic representation of the basin-wide SRM efficiency during calibration and validation period for the Hunza basin is presented in Figure 4b.The value of NS coefficient for the calibration and validation period was found to be 0.86 (Figure 4b).The measured and simulated flow curves indicated that SRM simulation at high and low flows was efficient during the whole study period.Pearson coefficient value was 0.92 for the correlation between the measured and simulated flow over the entire data period (Figure 5b) in the Hunza basin.It was also noticed that SRM overestimated the mid-level flows and underestimated the high flows in the Hunza basin (Figure 5b).(2000)(2001)(2002)(2003)(2004).The solid and dash lines represent regression and 1:1 lines respectively.

Simulation of the River Discharge by SRM Using the Zone-Wise Approach
In the zone-wise application, SRM was calibrated and validated over the same data period as in case of the basin-wide approach.The estimated values of parameters for the zone-wise simulation using SRM in the Astore and Hunza basins are presented in Table 6.

Simulation of the River Discharge by SRM Using the Zone-Wise Approach
In the zone-wise application, SRM was calibrated and validated over the same data period as in case of the basin-wide approach.The estimated values of parameters for the zone-wise simulation using SRM in the Astore and Hunza basins are presented in Table 6.
The efficiency of the SRM in terms of statistical analysis (using same statistical tests as in case of the basin-wide approach) with zone-wise application in the Astore and Hunza River basins over 5 years (2000)(2001)(2002)(2003)(2004) is presented in Table 7.For the Astore River basin, the application of the model was best in validation year 2002 for the zone-wise approach.The validation for year 2002 revealed NS coefficient value of 0.92 and volume difference of −9.8% (Table 7).The value of NS coefficient was not less than 0.80 which shows a satisfactory efficiency of SRM in the Astore basin.For Pearson and Kendall, rank correlation coefficient was not less than 0.91 and 0.75, respectively, in the Astore basin.
In Hunza River basin, the SRM model showed best efficiency for the calibration year 2000, NS coefficient was 0.90 with 0.4% volume difference (Table 7).NS value was never less than 0.81 which shows a good efficiency of SRM.For Pearson and Kendall rank correlation, the coefficient values were not less than 0.93 and 0.76, respectively, in the Hunza basin (Table 7).
Overall graphic representation of the zone-wise efficiency of the model during calibration and validation years for the Astore basin is presented in Figure 6a.The NS coefficient value for the calibration and validation period was found to be 0.87 and 0.85, respectively (Figure 6a).Measured and simulated flow curves suggest that the model simulated the high and low flows very efficiently during almost the entire data period except for the year 2003 where the simulated flow was largely under-estimated during early summer months.Pearson coefficient was 0.92 for the correlation between the measured and simulated flow over the entire data period (Figure 5c) in the Astore basin.Despite good overall zone-wise efficiency, SRM seems to be slightly over-estimating the high flows in the Astore basin (Figure 5c).Graphic representation of the zone-wise efficiency of the model during calibration and validation years for the Hunza basin is presented in Figure 6b.The NS coefficient value for the calibration and validation period was found to be 0.88 and 0.83, respectively (Figure 6b).Measured and simulated flow curves suggest that the SRM simulated the high and low flows very efficiently during almost the entire data period.The value of Pearson coefficient was 0.93 for the correlation between the measured and simulated flow over the entire data period (Figure 5d) in the Hunza basin.Scatter plot (Figure 5d) indicated a very good agreement between the low, mid and high level measured and simulated flow in the Hunza basin by applying zone-wise approach.and simulated flow curves suggest that the SRM simulated the high and low flows very efficiently during almost the entire data period.The value of Pearson coefficient was 0.93 for the correlation between the measured and simulated flow over the entire data period (Figure 5d) in the Hunza basin.Scatter plot (Figure 5d) indicated a very good agreement between the low, mid and high level measured and simulated flow in the Hunza basin by applying zone-wise approach.
Altogether the efficiency of the zone-wise simulation was greater as compared to the basin-wide simulation in the Hunza basin, but the later was superior in the Astore basin (Figure 5).The basinwide and zone-wise simulated discharges over the 2000-2004 period were plotted against each other to assess their conformity and it indicated a good correlation with a Pearson correlation coefficient value of 0.94 and 0.89 for for the Astore and Hunza basins, respectively (Figure 7).Altogether the efficiency of the zone-wise simulation was greater as compared to the basin-wide simulation in the Hunza basin, but the later was superior in the Astore basin (Figure 5).The basin-wide and zone-wise simulated discharges over the 2000-2004 period were plotted against each other to assess their conformity and it indicated a good correlation with a Pearson correlation coefficient value of 0.94 and 0.89 for for the Astore and Hunza basins, respectively (Figure 7).

Impact of Future Climate Scenarios on the Astore and Hunza River Runoff
The effect of climate change on the Astore and Hunza River runoff was estimated from the SRM by applying the RCP climate scenarios discussed in methodology part for mid-and late-21st century.The basin-wide simulated discharge by SRM in the year 2000 was used as the reference or base-year flow to project the future flows.
According to the SRM simulations under RCP 2.6, 4.5, and 8.5 for mid-21st century (2046-2065), there was approximately 13%, 20%, and 29% increase, respectively, in the mean River discharge in comparison with the base-year mean discharge (128 m 3 /s) in the Astore basin (Figure 8a).There was an approximate increase of 21%, 29%, and 44% in the mean discharge in comparison to base-year mean discharge (281 m 3 /s) under the RCP 2.6, 4.5 and 8.5, respectively, for mid-21st century (2046-2065) in the Hunza basin (Figure 8c).The results obtained after applying RCP 2.6, 4.5 and 8.5 for late-21st century showed an increase of 11%, 26% and 58% (Figure 8b) in the mean discharge, respectively, in the Astore basin.The results obtained after applying RCP 2.6, 4.5 and 8.5 for late-21st century showed an increase of 14%, 39% and 90% (Figure 8d) in the mean discharge, respectively, in the Hunza basin.Mean monthly discharges seem to experience an earlier shift of 15 days to one month in both the basins due to the earlier expected snowmelt under RCP 2.6, 4.5 and 8.5 for the mid-and late-21st century scenarios (Figure 8).
Considering only the summer high river flows, there was a significant increase in the mean summer discharge of the Astore and Hunza basin under RCP 2.6, 4.5 and 8.5 for the mid-as well as the late-21st century.A mean increase of ≈10-60% was projected in the summer flows of the Astore and Hunza rivers under RCP scenarios for mid-and late-21st century.This might be due to increased runoff from snow-and glacier-melt under increasing temperature as well as increasing rainfall trends as indicated in these scenarios.Similar trends in river discharge were projected until 2100 from Baltoro glacier melt under RCP 4.5 and RCP 8.5 scenarios [64].A clear warming trend in the UIB was indicated by [65] in all seasons as a result of global change in temperature and predicted rising river discharge due to glacial melting and increased rainfall runoff.But if these scenarios come true, the glaciers might vanish after generating excessive river flows up to a certain period.Similarly increasing discharge was projected under all RCP scenarios until the mid of this century and a slight decrease thereafter when ice cover (and ice thickness) will decrease considerably [66].Increasing temperature may lead to an accelerated ice melting which may increase the frequency and intensity of floods.Snowmelt predicted to occur earlier, with increasing ice thinning and even snow disappearing below 4000 m until 2100 [66].

Conclusions
Simulating stream flows from meltwater sustained sub-basins of the UIB are prerequisite for future water resource planning in Pakistan.The aim of this study was to simulate the runoff in the Astore and Hunza River basins under RCP climatic projections.The paucity of in-situ data in the difficult and remote terrain urged the use of MOD10A2 snow cover data with the accessible climatic data.Hydro-climatic datasets along with the topographical characteristics of the study area were used as input in the SRM to simulate the river flows.The following results were concluded: • SRM simulated the river flows efficiently during the calibration and validation years using the basin-wide and zone-wise approaches in the Astore and Hunza basins.Pearson correlation coefficient values of 0.95 and 0.92 were found for the basin-wide simulation in the Astore and Hunza basins, respectively, for the entire time period.Correlation values were 0.92 and 0.93 for the zone-wise simulations in the Astore and Hunza basins, respectively.• Application of SRM under climate change projections suggested an increase of 13-58% river flow for mid-to late-21 st century under different RCP scenarios in the Astore basin.This increase in river discharge was found to be 14-90% in case of Hunza basin.SRM was found to be an effective tool for simulation of snowmelt runoff in the data-scarce Karakoram and Western Himalayan regions targeted in this study.However, the availability of insitu precipitation data from high altitudes of the Astore and Hunza basins could improve the runoff simulations.The results indicated that if the climate change projections come true, an increasing discharge under all RCP scenarios could be generated until the end of this century and the seasoning of flows may also face a shift so that snowmelt starts earlier at least on the lower elevations.Increasing temperature may lead to accelerated snow and ice melting which may increase the frequency and intensity of floods.According to this study an efficient flood management plan may be needed; otherwise frequent floods might affect the already crippling economy.Therefore, assessing and mapping the flood-prone areas under future climatic projections may be useful to prevent future damages to human lives and infrastructures along with better planning and usage of water resources.

Conclusions
Simulating stream flows from meltwater sustained sub-basins of the UIB are prerequisite for future water resource planning in Pakistan.The aim of this study was to simulate the runoff in the Astore and Hunza River basins under RCP climatic projections.The paucity of in-situ data in the difficult and remote terrain urged the use of MOD10A2 snow cover data with the accessible climatic data.Hydro-climatic datasets along with the topographical characteristics of the study area were used as input in the SRM to simulate the river flows.The following results were concluded: • SRM simulated the river flows efficiently during the calibration and validation years using the basin-wide and zone-wise approaches in the Astore and Hunza basins.Pearson correlation coefficient values of 0.95 and 0.92 were found for the basin-wide simulation in the Astore and Hunza basins, respectively, for the entire time period.Correlation values were 0.92 and 0.93 for the zone-wise simulations in the Astore and Hunza basins, respectively.• Application of SRM under climate change projections suggested an increase of 13-58% river flow for mid-to late-21st century under different RCP scenarios in the Astore basin.This increase in river discharge was found to be 14-90% in case of Hunza basin.
SRM was found to be an effective tool for simulation of snowmelt runoff in the data-scarce Karakoram and Western Himalayan regions targeted in this study.However, the availability of in-situ precipitation data from high altitudes of the Astore and Hunza basins could improve the runoff simulations.The results indicated that if the climate change projections come true, an increasing discharge under all RCP scenarios could be generated until the end of this century and the seasoning of flows may also face a shift so that snowmelt starts earlier at least on the lower elevations.Increasing temperature may lead to accelerated snow and ice melting which may increase the frequency and intensity of floods.According to this study an efficient flood management plan may be needed; otherwise frequent floods might affect the already crippling economy.Therefore, assessing and mapping the flood-prone

Figure 1 .
Figure 1.Location map and digital elevation models of the study area (Hunza and Astore basins).

Figure 1 .
Figure 1.Location map and digital elevation models of the study area (Hunza and Astore basins).

Water 2019 , 19 Figure 2 .
Figure 2. Hypsometric curve of Astore and Hunza estimated through the digital elevation model.

Figure 2 .
Figure 2. Hypsometric curve of Astore and Hunza estimated through the digital elevation model.Mean monthly temperatures from 2000-2013 of Astore and Hunza basins are presented in Figure 3a.Mean monthly minimum temperature for the Astore (−1.9 • C) and Hunza (−11.7 • C) basins was observed in January while mean monthly maximum temperature was observed in July (11.5 • C) and August (20.7•C) for Hunza and Astore, respectively (Figure3a).Figure3bshows the monthly precipitation totals recorded in the Astore and Hunza River basins from 2000 to 2013.The Figure3bshows that maximum precipitation was received during the September-April, in winter and spring season, as snow, because the temperature during these months was below or near freezing (Figure3a) in both basins.Most of this precipitation originates from Westerlies.Astore and Hunza rivers are gauged at Doyian and Dainyor hydrometric stations (Figure1), respectively.Astore river joins the Indus just downstream of its flow gauge.The Hunza river joins the Indus River few kilometers downstream after the confluence with the Gilgit River at Dainyor.River discharges converted to equivalent water depths (in mm) upon dividing by the basin area for both basins are presented in Figure3c, which represents that Astore has more runoff than the Hunza which is also indicative of the more precipitation amounts (Figure3b) in the Astore than Hunza basin.Astore also has higher temperatures than the Hunza (Figure3a), so it may have more melt runoff generation than the Hunza basin.Astore River flow starts increasing in the month of March, and a peak discharge is achieved in June, as observed over the available data period (Figure3c) when monsoon rainfall is superimposed by the peak melting season.The Hunza River flow starts increasing in the month of April and reaches its maximum during July as observed over the available data period (Figure3c).The reason for delayed flow rising limb in the Hunza River (in April as compared to March in Astore river) is the mean temperature which only starts to reach above freezing in the month of April (Figure3a).Mean snow cover area estimated using the MODIS snow product for the period of 2000-2016 is presented in Figure3d.Hunza basin has slightly higher mean annual SCA (64%) than the Astore basin (57%).Maximum mean snow cover area is found to be in the March (end of snow accumulation season) and minimum in August (end of snowmelt season) in these basins (Figure3d).

Figure 3 .
Figure 3. Hydro-climatic data representation in the Astore and Hunza basins (a) Mean monthly temperature variation (average of 2000-2013), (b) Total monthly precipitation (average of 2000-2013 totals), (c) Mean monthly runoff (average of 2000-2013), and (d) Basin-wide mean monthly snow cover area (with standard deviation bars) over a period of 2000-2016.

Figure 3 .
Figure 3. Hydro-climatic data representation in the Astore and Hunza basins (a) Mean monthly temperature variation (average of 2000-2013), (b) Total monthly precipitation (average of 2000-2013 totals), (c) Mean monthly runoff (average of 2000-2013), and (d) Basin-wide mean monthly snow cover area (with standard deviation bars) over a period of 2000-2016.

Figure 5 .
Figure 5. Correlation between observed and simulated discharge (Q) is presented in (a) and (b) panels for the basin-wide application and (c) and (d) panels for the zone-wise application of SRM in the Astore and Hunza basins over study period(2000)(2001)(2002)(2003)(2004).The solid and dash lines represent regression and 1:1 lines respectively.

Figure 5 .
Figure 5. Correlation between observed and simulated discharge (Q) is presented in (a) and (b) panels for the basin-wide application and (c) and (d) panels for the zone-wise application of SRM in the Astore and Hunza basins over study period(2000)(2001)(2002)(2003)(2004).The solid and dash lines represent regression and 1:1 lines respectively.

Figure 7 .
Figure 7. Agreement between the basin-wide and zone-wise SRM simulation in the (a) Astore and (b) Hunza basins.The solid and dash lines represent regression and 1:1 lines respectively.

Table 1 .
Parameters and variables used in the Snowmelt-Runoff Model (SRM).

Table 2 .
Representative Concentration Pathways (RCPs) scenarios downscaled for the Indus River Basin and used in the current study to project future river flows.

Table 3 .
Characteristics of elevation zones estimated from DEM for Astore and Hunza basin for the zone-wise application of SRM.

Table 4 .
Parametric values used in basin-wide application of SRM for the study area.

Table 5 .
Efficiency of the basin-wide SRM application to simulate the Astore and Hunza river discharge.

Table 6 .
Parametric values used in the zone-wise application of the SRM for the Astore and Hunza basins.

Table 6 .
Parametric values used in the zone-wise application of the SRM for the Astore and Hunza basins.

Table 7 .
Efficiency of the zone-wise SRM application to simulate Astore and Hunza river discharge.