Di ﬀ erentiating Snow and Glacier Melt Contribution to Runo ﬀ in the Gilgit River Basin via Degree-Day Modelling Approach

: In contrast to widespread glacier retreat evidenced globally, glaciers in the Karakoram region have exhibited positive mass balances and general glacier stability over the past decade. Snow and glacier meltwater from the Karakoram and the western Himalayas, which supplies the Indus River Basin, provide an essential source of water to more than 215 million people, either directly, as potable water, or indirectly, through hydroelectric generation and irrigation for crops. This study focuses on water resources in the Upper Indus Basin (UIB) which combines the ranges of the Hindukush, Karakoram and Himalaya (HKH). Speciﬁcally, we focus on the Gilgit River Basin (GRB) to inform more sustainable water use policy at the sub-basin scale. We employ two degree-day approaches, the Spatial Processes in Hydrology (SPHY) and Snowmelt Runo ﬀ Model (SRM), to simulate runo ﬀ in the GRB during 2001–2012. The performance of SRM was poor during July and August, the period when glacier melt contribution typically dominates runo ﬀ . Consequently, SPHY outperformed SRM, likely attributable to SPHY’s ability to discriminate between glacier, snow, and rainfall contributions to runo ﬀ during the ablation period. The average simulated runo ﬀ revealed the prevalent snowmelt contribution as 62%, followed by the glacier melt 28% and rainfall 10% in GRB. We also assessed the potential impact of climate change on future water resources, based on two Representative Concentration Pathways (RCP) (RCP 4.5 and RCP 8.5). We estimate that summer ﬂows are projected to increase by between 5.6% and 19.8% due to increased temperatures of between 0.7 and 2.6 ◦ C over the period 2039–2070. If realized, increased summer ﬂows in the region could prove beneﬁcial for a range of sectors, but only over the short to medium term and if not associated with extreme events. Long-term projections indicate declining water resources in the region in terms of snow and glacier melt.


Study Region Description
The Gilgit basin is a sub-basin of UIB and encompasses the region from 35.8-37 • N to 72.5-74.4 • E, including the eastern part of the Hindukush range. The drainage area of this basin incorporates 12,000 km 2 with an elevation range of 1481-7134 m.a.s.l, as shown in Figure 1. The primary controls on precipitation in the basin are associated with disturbances in the upper-level westerlies and the summer monsoon [13]. The basin has approximately 685 glaciers and 605 glacial lakes. Eight of these lakes have been identified as a potential threat due to Glacial Lake Outburst Floods (GLOF) [27]. The careful selection of GRB involves its cryospheric extent based on snow and glacier coverage. Dominant snow coverage makes this basin perfect for using SRM, while slight glacier coverage in terms of clean-ice and debris-covered glaciers ensure smooth application of SPHY.

Study Region Description
The Gilgit basin is a sub-basin of UIB and encompasses the region from 35.8-37° N to 72.5-74.4° E, including the eastern part of the Hindukush range. The drainage area of this basin incorporates 12,000 km 2 with an elevation range of 1481-7134 m.a.s.l, as shown in Figure 1. The primary controls on precipitation in the basin are associated with disturbances in the upper-level westerlies and the summer monsoon [13]. The basin has approximately 685 glaciers and 605 glacial lakes. Eight of these lakes have been identified as a potential threat due to Glacial Lake Outburst Floods (GLOF) [27]. The careful selection of GRB involves its cryospheric extent based on snow and glacier coverage. Dominant snow coverage makes this basin perfect for using SRM, while slight glacier coverage in terms of clean-ice and debris-covered glaciers ensure smooth application of SPHY.

Hydro-Climatological Data
Three organizations, the Pakistan Meteorological Department (PMD), the Surface Water Hydrology Project Water and Power Development Authority (SWHP-WAPDA), and the Snow and Ice Hydrology Project (SIHP-WAPDA), are responsible for the network of meteorological and hydrological stations installed within Pakistan's territory [38]. World Meteorological Organization (WMO) guidelines are followed for the installation of instrumentation, data acquisition, and subsequent distribution of data to users [39]. The PMD operates six manual valley-based stations, including Gilgit (1460 m.a.s.l) and Gupis (2156 m.a.s.l) and data for these stations are available from the 1960s onwards. The streamflow network is primarily managed by SWHP-WAPDA. High altitude data collection is operated by SIHP-WAPDA, which consists of 12 high altitude Automatic Weather Stations (AWS) including Ushkore (3051 m.a.s.l) and Yasin (3280 m.a.s.l) [13,40]. Daily hydroclimatic data (e.g., precipitation, temperature, and flows) at all four stations, including a streamflow gauge on Gilgit River at Gilgit (Figure 1), were provided by the PMD and WAPDA.

Meteorological and Climate Model Data
The SRM and SPHY models use daily temperature and precipitation data to simulate runoff. The required meteorological data are available from Future Water, for the Hindukush-Himalayan Region, as shown in Table 1. These meteorological data are derived from the Watch Forcing Data ERA-INTERIM (WFDEI). We subsequently converted the available daily Tmax, Tmin, Tavg, and precipitation (Prec) data into the required format for input into SPHY and SRM.

Hydro-Climatological Data
Three organizations, the Pakistan Meteorological Department (PMD), the Surface Water Hydrology Project Water and Power Development Authority (SWHP-WAPDA), and the Snow and Ice Hydrology Project (SIHP-WAPDA), are responsible for the network of meteorological and hydrological stations installed within Pakistan's territory [38]. World Meteorological Organization (WMO) guidelines are followed for the installation of instrumentation, data acquisition, and subsequent distribution of data to users [39]. The PMD operates six manual valley-based stations, including Gilgit (1460 m.a.s.l) and Gupis (2156 m.a.s.l) and data for these stations are available from the 1960s onwards. The streamflow network is primarily managed by SWHP-WAPDA. High altitude data collection is operated by SIHP-WAPDA, which consists of 12 high altitude Automatic Weather Stations (AWS) including Ushkore (3051 m.a.s.l) and Yasin (3280 m.a.s.l) [13,40]. Daily hydroclimatic data (e.g., precipitation, temperature, and flows) at all four stations, including a streamflow gauge on Gilgit River at Gilgit (Figure 1), were provided by the PMD and WAPDA.

Meteorological and Climate Model Data
The SRM and SPHY models use daily temperature and precipitation data to simulate runoff. The required meteorological data are available from Future Water, for the Hindukush-Himalayan Region, as shown in Table 1. These meteorological data are derived from the Watch Forcing Data ERA-INTERIM (WFDEI). We subsequently converted the available daily T max , T min , T avg , and precipitation (Prec) data into the required format for input into SPHY and SRM. We employed global daily downscaled data (GDDP) for the assessment of climate flows due to a changing climate. This dataset is provided for assessing climate change impacts at local to regional scales and to enhance public understanding of possible future climate patterns at the spatial scale of individual towns, cities, and watersheds. The GDDP dataset contains downscaled climate scenarios for the globe derived from General Circulation Model (GCM) runs conducted as part of the Coupled Model Intercomparison Project Phase 5 (CMIP5) [41]. These datasets provide two of four greenhouse gas emissions scenarios, referred to as Representative Concentration Pathways (RCPs), namely RCP 4.5 and RCP 8.5 [42]. The CMIP5 GCM simulations were undertaken in support of the Intergovernmental Panel on Climate Change's (IPCC) Fifth Assessment Report (FAR) [43]. A number of limitations exist with the use of these data. Firstly, regional or local climate analyses typically require detailed spatial and temporal climate information, and consequently, the coarse grid resolution may act to smooth some of the local climate variability. Secondly, global projections, particularly at the regional scale, display (systematic) biases which need to be accounted for before use [44].
While recognizing these shortcomings, we employed climate data from the Can-ESM2 and Nor-ESM1-M Global Climate Models (GCMs), based on RCP 4.5 and 8.5, for use in SRM within GRB, to assess the potential impact of climate change on future river flows in the region. The NEX-GDDP datasets are generated by the Bias-Correction Spatial Disaggregation (BCSD) method featuring a statistical downscaling approach to address limitations of directly employing GCM output [45][46][47]. The bias correction method employs spatial details acquired from observational datasets for the interpolation of GCM outputs to higher resolution grids. The algorithms compare GCM output with observed climate data for a common or overlap period; the difference is then used to statistically Atmosphere 2020, 11, 1023 6 of 25 'bias correct' or adjust the climate model projections based on the assumption that any differences are systematic. The downscaled datasets generated by BCSD also assume that the relative spatial patterns in precipitation and temperature during the historical period will remain constant under the future climate change scenarios-referred to as pattern scaling. The relative spatial patterns are then used to disaggregate the bias-corrected climate model data to the same grid resolution as the observations, using a bilinear interpolation technique. The spatial details of fine-grained observations are preserved using this technique, for example, by computing observed daily precipitation values which are multiplied by a scaling factor to adjust the daily interpolated GCM projections to produce pattern scaled, regional/local climate projections [48].

Ancillary/Satellite Derived Data
SRM requires snow cover products as a primary input; due to the scarcity of available in situ data, remotely sensed data are typically utilized [20]. MODIS 8 days L3 Global 500 m grid processed snow cover products (MOD10A2) were used to calculate the extent of snow cover area during 2000-2018 within the GRB. We processed 430 available images, from April 2000 to September 2018, using the maximum extent of snow cover spanning the 8-day coverage to minimalize error due to cloud cover. The selected images were subsequently mosaicked and projected using ENVI + IDL 4.7; the images for the Gilgit basin were then extracted from these mosaicked scenes. Finally, we calculated the snow cover area (SCA) for seven different altitudinal zones within GRB. However, the lower albedo prevents the MODIS platform from acquiring precise debris-covered glacier fractions.
Daily forcing data are available from the SPHY database for the Hindukush-Himalayan region; alternatively, user-defined data can be employed. In this study, we employed a higher resolution digital elevation model (DEM) derived from NASA's Shuttle Radar Topography Mission (SRTM) obtained at 90 m resolution. All user-defined datasets were subsequently converted to comma delimited (e.g., CSV) format using the method outlined in the SPHY manual. We used similar CSV formatted files for the basic meteorological forcing for the stations outlined in Table 1.

Snowmelt Runoff Model (SRM)
SRM is a conceptual, deterministic degree-day model used to simulate daily streamflow values resulting from snowmelt and precipitation in mountainous regions. The basic inputs for the SRM are daily temperature, precipitation, and SCA.
Streamflow is calculated as Q n+1 = [C sn · a n + (T n + ∆T n ) C Rn P n ] A.1000 86400 where Q n+1 is the discharge at day n+1 (m 3 /s); C sn is the snow runoff coefficient; a n is the degree-day factor on day n (cm/ • C/day); T n + ∆T n are the degree days ( • C); S n is the fractional snow cover; C Rn is the rain runoff factor; P n is the rain on day n (cm); A is total area(km 2 ); K n+1 is discharge recession coefficient. An illustration of the methodology is shown in Figure 2.    A detailed description of the SRM model is provided by [20]. The SRM is more sensitive to snow cover data as compared to precipitation data. Therefore, it is preferable to use snow cover data, when available precipitation data are sparse or less reliable [24]. Furthermore, SRM is unable to calculate glacier area coverage by default; an assumption is made that the glacier extent is persistent for the calculation of actual contribution from glaciers for debris-covered as well as clean-ice glaciers. This is based on the fact that glacier area changes occur over longer durations as compared to the seasonal snow cover which melts seasonally. Another limitation of the model is introduced at the basin scale, where the use of a single DDF value is assumed to be representative of the entire basin and for both debris-covered and clean-ice glacier fractions. An alternative approach involves the spatial disaggregation of the basin into representative zones (zone-wise application), depending upon the basin size [49]. Here, we used SRM in zone-wise application by extracting seven altitudinal zones from the entire GRB ( Figure 3). We divided the basin zones using a DEM and subsequently calculated the SCA for each zone. The study area of GRB calculated by SPHY was around 12,300 km 2 which is approximately equal to the area calculated using the DEM. The zones are outlined in Table 2.   An overview of the parameters and calibrated values for SRM is provided in Table 3. TLR is considered as the most sensitive parameter in the temperature-based/degree-day modelling approach. Previous studies [50] reported that underestimation of the TLR may lead to an underestimation of melting rates and vice versa. TLR typically varies between −0.98 to 0 • C/100 m, referred to as dry adiabatic and isothermal lapse rate, respectively, and the global standard lapse rate of −0.65 • C/100 m is widely used [51]. Previous studies within high altitude areas [24,52,53] have employed annual TLR values ranging from −0.2 to −0.9 • C/100 m in model simulations. We used station data up to 3280 m (Yasin) as the temperature input for model runs. We extrapolated the temperature values above this elevation range in different altitudinal zones using a lapse rate ranging from −0.42 to −0.59 • C/100 m for various zones. Previously, [24,28] used −0.48 to −0.76 and −0.64 • C/100 within the Hunza and Gilgit basins, respectively, which is slightly higher than the values used in adjacent regions such as the Tibetan Plateau (TP) [26]. The authors also proposed that the constant value of TLR does not accord with the environmental conditions experienced in the TP and adjoining regions. Consequently, we employed a TLR value ranges −0.42 to −0.53 • C/100 m, suggested by [54], which is closer to the average global TLR. Table 3. Calibration parameter values for zone-wise SRM application to the Gilgit River basin.

Parameters
Parameters for Each Altitudinal Zones 1(1460-2500) 2(2500-3500) 3(3500-4000) 4(4000-4500) 5(4500-5000) 6(5000-5500) 7(5500-7150) The melting rates for clean ice and snow vary over UIB. According to [55], the melt rate is 0.7 and 0.4 cm/ • C/day for clean ice and snow, respectively. Previous studies reported differences in DDF with regards to elevation, while [47] reported a depressed DDF attributed to a lower value associated with the density of fresh snow. Some authors [24] used 0.5 to 0.7 cm/ • C/day for zone-wise application, however, [28] used 0.1 to 0.6 cm/ • C/day during early to late ablation period. Some authors [26] have employed a constant DDF value for snow and glacier ice for the TP. They also reported that variations in DDF result in an overestimation of new snow, resulting in faster snowmelt in the model simulation. These authors suggested the use of a varying DDF based on season and elevation zones. Here, we employ a varying DDF for each month, from April to September i.e., 0.5 for April-May, 0.75 for June-August and 0.35 cm/ • C/d for September. These values were calibrated based on geodetic mass balance data for snow, clean-ice and debris-covered glaciers in the Karakoram region [37]. The GRB is a snow dominant basin, therefore DDF is lower in the beginning but increases during peak ablation period for melting of the snow and the clean-ice glacier. However, it decreases in September when melt contribution generates by the debris-covered glacier. This value can be opted out depending upon the glacier extent (clean-ice/debris-covered) for rest of the basins in UIB. We also used the varying value for Rainfall contributed area (RCA) as '1 for June-August and '0 for May and September. The duration of the period concerned was April to September i.e., snowmelt/glacier melt period. It is assumed that when rain falls on previous snow cover, an equal amount of water will be released from the snowpack so that the rain from that area could be added as snowmelt [56]. The melting effect of the rain was neglected due to the small amount of additional heat supplied by the liquid precipitation [57], ultimately a pragmatic simplification. Some other parameters such as rainfall (Cs) and snow (Cr) runoff coefficients which dictate the conversion of rainfall and snowfall into the runoff. Similarly, recession coefficients (Xc and Yc) account for the runoff contribution by snowmelt occurred on the previous day. These parameters were calibrated according to previous studies as suggested by [28]. Calibration parameter values for the zone-wise application of SRM are given in Table 3.

Spatial Processes in Hydrology (SPHY)
The SPHY is a spatially distributed, leaky bucket model which is applied in cell by cell mode [35]. The model uses high-resolution cryospheric-hydrological characteristics and is available as an open-source versions (2.0-2.2) for users [37]. The model is written in Python using PCRaster as the dynamic modelling framework [58]. The SPHY was developed by Future Water [58] for runoff simulation for both the present and future, assuming available data. SPHY can be used for understanding hydrological processes such as glacier processes [36].
It is a user-friendly model which requires more extensive inputs in addition to daily temperature and precipitation data. The model is more complex, theoretically and practically, in terms of data inputs and outputs compared to SRM. The most important step in the SPHY simulation involves the selection of crop coefficient (Kc) values for the study area which can be estimated according to the land use classes. However, Kc values are known to vary over time and space [36]; therefore, it is preferable to use a variable (e.g., seasonal; spatial) crop coefficient rather than using a single static value. The transportation of water through the delineated river network is affected by recession coefficient (Kx) [36], which must be opted out between 0 and 1. Here, we employ the same calibrated DDF based on geodetic mass balance data [37] for snow, clean-ice and debris-covered glaciers. The duration of percolated water within the second to the third layer to achieve the groundwater store is controlled by delta Groundwater (deltaGw), while alpha Groundwater (alphaGw) controls baseflow recession curve. The values of delta Gw can be opted out from 1 to several days, however, Gw ranges between 0 to 1 [59]. During an initial run, SPHY generates land use information which can be readily incorporated for subsequent processing. Parameters values used in SPHY are shown in Table 4. SPHY was run on a daily time step and spatial resolution of 500 m, covering the entire basin. The basic methodology of runoff routing and melt modules employed in SPHY are outlined below.
SPHY estimates total runoff using the following relationship where QTot denotes total runoff in (mm) for a particular cell, RRo represents rainfall runoff in (mm), SRo stands for snow runoff (mm), GRo is glacier runoff (mm), and BF is baseflow (mm) which can be from the third soil layer or lateral flow from the second soil layer. The glacier melt is calculated by the degree day modelling approach [60]. The glaciers respond differently depending upon whether they are debris-covered or have clean-ice due to the supraglacial debris layer modifying the interaction between the glacier surface and the overlying atmosphere [61]. The SPHY uses different values of the degree day factor (DDF) for clean-ice and debris-covered glaciers. The equation is where ACI stands for daily melt rate (mm) for clean-ice glacier, DDF CI is a calibrated degree day factor (mm/ • C/day), and F CI represents fraction of clean-ice glacier within fractional glacier cover (GlacF). The GlacF is calculated when the model creates module of SPHY. The melting rate from a debris-covered glacier is calculated by a similar equation, where ADC is the daily melt rate (mm) for a debris-covered glacier, DDFDC is a calibrated degree day factor (mm/ • C/day), and FDC represents the fraction of clean-ice glacier within fractional glacier cover (GlacF). The GlacF for a debris-covered glacier is calculated during the area selection step. The melt rate from debris-covered and clean-ice glaciers is summed up to obtain total melting rate, which is calculated by the equation.
The daily potential of snowmelt is calculated as follows where A pot,t is potential snowmelt (mm) on day t; DDFS is calibrated degree day factor for snow (mm/ • C/day). Actual snowmelt is calculated by the equation where A act,t is actual snowmelt (mm) on day t, which is limited by the snow storage at the end of the previous day. SS t−1 is snow storage on day t − 1. Snow storage form day t − 1 is updated to day t using actual snowmelt and solid precipitation. The fraction of snowmelt freezes within snowpack, which does not become a part of runoff immediately. It is added to the snow storage when the temperature falls below melting point on day t − 1 as follows where SS t (mm), SS t−1 , and SSW t−1 denote snow storage on day t, snow storage on day t − 1, and frozen meltwater on day t − 1, respectively; P s,t is solid snow on day t; A act, t is actual snowmelt on day t. When the temperature is above the melting point and there is no water to be frozen, snow runoff is generated, which is calculated as follows where SRo t is snow runoff in mm, ∆SSW is change in meltwater stored in snowpack and calculated as

Model Evaluation
Model accuracy was assessed using the coefficient of determination (R 2 ), given as The Q i denotes the daily observed flow, Q i ' is the daily simulated flow, and Q is the average constant stream of the data. The value of R 2 will be 1 when the root mean square prediction error approaches 0.
Similarly, the volume difference between observed and measured flow is calculated by the relation given as where V is measured and V' is the computed runoff volume. This relationship is used for accuracy calculation (WMO). An over or underestimation of flow is reflected in negative and positive values of Dv, respectively. We evaluated both models based on Nash-Sutcliffe Efficiency (NSE), R 2 , Root-Mean-Square-Error (RMSE) and bias values in achieving Goodness of Fit (GOF) during the calibration process [59].

Mann Kendall and Sen's Estimator
We used the Mann-Kendall [62] and Sen's slope estimator [63] for the magnitude and direction of climatic trends (precipitation, temperature, and snow cover) during the selected study period. These methods for trend identification have been widely use and established by a significant number of previous studies carried out within UIB, such as [13,38,39,64,65].

Results
Time series of precipitation and temperature from GRB, for the period from 2000 to 2014, are shown in Figure 4. The highest maximum monthly temperatures were recorded at Gilgit (Figure 4a

Results
Time series of precipitation and temperature from GRB, for the period from 2000 to 2014, are shown in Figure 4. The highest maximum monthly temperatures were recorded at Gilgit (Figure 4a), while the lowest minimum monthly temperatures were recorded at Yasin (Figure 4b). The maximum annual (Figure 4d) and winter precipitation (Figure 4f) observed at the Yasin and Ushkore stations were typically higher than those of Gupis and Gilgit. The maximum observed winter ( Figure 4f) and summer (Figure 4e) precipitation, of 450 and 348 mm, respectively, were recorded at Ushkore and Yasin during 2004 and 2010. We estimated the snow and glacier extent of GRB during the selected study period. According to this analysis, approximately 9% of the Gilgit basin was covered by glaciers. This finding is consistent with that of [6,27], who reported approximately 8% glacier coverage. The extent of cleanice, debris-covered ice, and annual snow persistency in percentage is illustrated in Figure 5. The clean-ice covered area is calculated as 945 km 2 , while the area of debris-covered glaciers is 145 km 2 , again consistent with [27]. We estimated the snow and glacier extent of GRB during the selected study period. According to this analysis, approximately 9% of the Gilgit basin was covered by glaciers. This finding is consistent with that of [6,27], who reported approximately 8% glacier coverage. The extent of clean-ice, debris-covered ice, and annual snow persistency in percentage is illustrated in Figure 5. The clean-ice covered area is calculated as 945 km 2 , while the area of debris-covered glaciers is 145 km 2 , again consistent with [27]. The gradual melting of snow cover during the melt season of 2012 is shown in Figure 6. During 2012, snowmelt begins in late spring (April) and continues throughout the early-summer months; the seasonal snow cover melt extent reaches maximum by the end of July. During August, glacier ice melt commences, except over high-altitude zones due to lower temperatures, and continues until late The gradual melting of snow cover during the melt season of 2012 is shown in Figure 6. During 2012, snowmelt begins in late spring (April) and continues throughout the early-summer months; the seasonal snow cover melt extent reaches maximum by the end of July. During August, glacier ice melt commences, except over high-altitude zones due to lower temperatures, and continues until late September. Only snow at high altitude persists throughout the ablation season. Accumulation initiates in the middle of October and continues up to March of the following year. Figure 7 illustrates the monthly depletion curve values for the period 2001-2018 in the Gilgit basin; we observed stable or slightly increasing (significant) Cumulative Depletion Curve (CDC) values during the 2001-2018 period, based on the MK test. We also observed an insignificant positive trend (less than 90% confidence interval) in the SCA of GRB during the months of July to September.
Atmosphere 2020, 11, x FOR PEER REVIEW 14 of 24 subsequently decreased during 2011 and 2012 due to increased contributions from the snow and glacier meltwater. We also noted the changes in SCA during the ablation period in all elevation zones as illustrated in Figure 9b. It is evident that the snow remains intact (low variation) in high elevation zones at the beginning of ablation period and melting appears in the form of slight variations during the peak ablation period (Figure 9c). However, melting appears stronger in the remaining zones as CDC reduces to 70% for Zone 6 during peak ablation. The highest reduction, exhibited in Zone 3, resulted in the complete melt out of the snow by the end of the ablation period. The CDC values are based on the maximum snow cover extent acquired weekly from April to September. We interpolated these values to derive daily snow extent values; the higher CDC values reflect the maximum extent of SCA. We observed that most of the changes appeared in the lower elevation zones in the form of rapid melting, whereas the high elevation zones featured reduced melting.     The SRM exhibits daily snowmelt and rainfall values (cm) in each zone as output. The snowmelt in each zone represents both snow and glacier melt contributing by each zone. Due to lower accumulation of snow and permanent glacier ice, the melting rate in the lower elevation zones remains slow, so all of the snow accumulated in this zone melts by the end of May while melting over higher elevation zones (5, 6, and 7) continues up to September, as shown in Figure 8. The flows during the rest of the year (autumn, winter, and spring) remain low compared to summer. The SRM determines whether precipitation is treated as liquid or solid, depending on the value of T crit . We noted that a portion of the precipitation is consequently attributed to snowmelt depth, while the remainder contributes to flow.   The contribution of rainfall in the lower elevation zones is high and gradually decreases until it reaches a negligible amount in the upper regions of the high-altitude zones. Similarly, the snow and glacier melt contribution remains low in the lower elevation zones, particularly during peak ablation season. The increasing dominance of the snow and glacier melt water contribution with increasing elevation is also marked in the average values calculated over the 2001-2012 period, as shown in Figure 9a. The average rainfall contribution in Zone 1 is 80% and decreases to 8% in Zone 7. Similarly, the snow and glacier contribution is 18% in Zone 1 and increases to a maximum level of 92% in the uppermost zone, Zone 7. We also observed that the maximum rainfall contribution to flow occurred during 2008 and 2010, reaching 37% in 2010. This percentage contribution is consistent with the increase in precipitation observed at all stations in GRB during 2010. The rainfall contribution subsequently decreased during 2011 and 2012 due to increased contributions from the snow and glacier meltwater. We also noted the changes in SCA during the ablation period in all elevation zones as illustrated in Figure 9b. It is evident that the snow remains intact (low variation) in high elevation zones at the beginning of ablation period and melting appears in the form of slight variations during the peak ablation period (Figure 9c). However, melting appears stronger in the remaining zones as CDC reduces to 70% for Zone 6 during peak ablation. The highest reduction, exhibited in Zone 3, resulted in the complete melt out of the snow by the end of the ablation period. The CDC values are based on the maximum snow cover extent acquired weekly from April to September. We interpolated these values to derive daily snow extent values; the higher CDC values reflect the maximum extent of SCA. We observed that most of the changes appeared in the lower elevation zones in the form of rapid melting, whereas the high elevation zones featured reduced melting.    The flow contribution using SPHY is illustrated in Figures 10 and 11. We found that SPHY is more capable in terms of NSE and R 2 in order to accomplish GOF during calibration, relative to SRM, with the exception of the values for RMSE and bias which are slightly higher for both models (Figure 10a). Similarly, the SPHY can determine the glacier contribution as a distinct component of the total melt (Figure 10b), while SRM indicates only rain and combined snow and glacier contribution. The average simulated runoff revealed the prevalent snowmelt contribution as 62%, followed by the glacier melt 28% and rainfall 10% in GRB. We also found an increasing rainfall contribution over the period 2001-2012 in both models. Similarly, the snow and glacier contribution assessed by SRM and SPHY is more than 60% during the same period, with the exception of the years 2005 and 2006, where SRM shows a greater contribution from rainfall. While the snowmelt contribution remains stable, the glacier contribution was found to increase slightly over the 2001-2012 period, as shown in Figure 10b. Within SPHY, rainfall runoff is calculated on the basis of gridded precipitation data at high altitude and is generated in the form of direct/immediate runoff at the hill slope (Figure 11a). Snow runoff is calculated in the form of delayed runoff contributing after snowmelt has occurred, during the ablation season at Gilgit at Gilgit River, as shown in Figure 11b. Compared to snow runoff, the glacier contribution is lower, based on the glacier extent and the contributed runoff (Figure 11c). glacier contribution was found to increase slightly over the 2001-2012 period, as shown in Figure 10b. Within SPHY, rainfall runoff is calculated on the basis of gridded precipitation data at high altitude and is generated in the form of direct/immediate runoff at the hill slope (Figure 11a). Snow runoff is calculated in the form of delayed runoff contributing after snowmelt has occurred, during the ablation season at Gilgit at Gilgit River, as shown in Figure 11b. Compared to snow runoff, the glacier contribution is lower, based on the glacier extent and the contributed runoff (Figure 11c).   glacier contribution was found to increase slightly over the 2001-2012 period, as shown in Figure 10b. Within SPHY, rainfall runoff is calculated on the basis of gridded precipitation data at high altitude and is generated in the form of direct/immediate runoff at the hill slope (Figure 11a). Snow runoff is calculated in the form of delayed runoff contributing after snowmelt has occurred, during the ablation season at Gilgit at Gilgit River, as shown in Figure 11b. Compared to snow runoff, the glacier contribution is lower, based on the glacier extent and the contributed runoff (Figure 11c).   We subsequently analyzed the potential impact of climate change on river flows using SRM based on two RCPs (4.5 and 8.5) employed by the Can-ESM2 and Nor-ESM1-M GCMs, for the entire 21st century within GRB. These two GCMs were selected on the basis of previously published research on this region [28,36]. The future period was divided into three data series, where the first series overlapped the present period as T 1 (2010-2039), T 2 (2040-2069), and T 3 (2070-2100). The reference years of 2010, 2011, and 2012 were used as the baseline to assess future flows. The results are outlined in Table 5. Trend is significant at 99.99% (***), 99% (**) and 95% (*) confidence interval. Confidence interval ≤90% = Insignificant trend.
Based on RCP 4.5 from the Can-ESM2 and Nor-ESM1-M models, increases in average annual temperature of 1.7 and 0.7 • C are projected, respectively, resulting in an increase in flows of 10.87% and 5.6% between T 1 and T 3 . However, an increase of 1.3 • C increased the annual flow by 10.6% during T 2 . Similarly, for RCP 8.5, temperatures are projected to increase by 2.4, 2.1, and 2.6 • C resulting in increased flows at the rates of 16.4% 14.5%, and 19.08%, respectively, during T 1 , T 2 , and T 3 . We observed the similar increases in summer temperature at the rates of 0.8 and 0.8 • C, resulting in enhanced flows at the rates of 5.81% and 6.1% during T 1 and T 3 , respectively, for RCP 4.5 and Nor-ESM1-M. Similarly, we found a positive trend of average summer temperature at the rate of 0.2 • C, causing an increased flow at 2.6%, but was statistically insignificant. We noted positive trends at the rates of 2.9, 1.9, and 1.8 • C causing increased flows at the rates of 18.2%, 11.3%, and 10.43%, respectively, during T 1 , T 2 , and T 3 .
We found slightly increased flows during T 1 for both scenarios, although an increase in flows was estimated particularly for RCP4.5 during T 2 ; similarly, we calculated the future projections of mean annual temperature • C per period from Can-ESM2 for summer in the Gilgit basin, as shown in Table 5. We found significant positive trends in temperature during T 1 , while annual and seasonal trends were revealed as insignificant during T 3 . Our analysis suggested significant positive trends during T 3 for RCP8.5, but these trends were exhibited as insignificant during T 1 and T 2 . We estimate that the average Gilgit River summer flows are likely to increase by 5.6% and 19.8%, associated with increases of 0.77 and 2.60 • C in the average summer temperatures during T 2 under (RCP) 4.5 and 8.5 scenarios, respectively.

Discussion
The present study enabled the assessment of snow cover trend and runoff simulations in addition to the differentiation of snow/glacier melt contribution to runoff within GRB. Some recent studies [34] suggested an increasing SCA trend within the Gilgit Basin attributed to decreased liquid precipitation with increased solid precipitation. Previous studies concluded a reduced SCA during winter and autumn but growing summer SCA within the Hunza [24] and Gilgit basins [6] during 2001-2012. The growing summer SCA is associated with winter warming and summer cooling effects [14]. Our findings of increasing SCA confirm that the increasing trend now extends to cover the period 2001-2018. The expansion of SCA is likely attributable to increasing precipitation trends [65] and significantly reduced temperature [38] at high altitude stations in GRB. Snow cover and snow accumulation amount are linked directly [6]; they concluded that the new snow remains for a shorter period at the steeper angle at high altitude. According to our findings, increasing snow cover is attributed to increased winter precipitation within this basin during the last two decades [39,66].
The glaciers in the GRB gain less height as compared to the Hunza River Basin, so the persistence of initial snow increases storage capacity to hold new (incoming) snow until melting temperature is achieved to initiate melting. We observed the minimum variation in SCA particularly in the highest elevation zone (5501-7150 m). These findings are also consistent with significant decreases in maximum temperature within the Gilgit basin at Gilgit, Gupis, Ushkore, and Yasin [38]. The presence of increasing cloud cover and increasing number of wet days [66] for a similar region is also in agreement with such increasing/stable snow cover. The consistent increase in cloud cover and frequent monsoonal incursions cause reduced downward radiation, ultimately resulting in a cooling environment [13]. The decreasing temperature in summer decreased the melting rate for new snow as well as for permanent snow and glaciers. The recent temperature and precipitation trends in the Gilgit basin were comprehensively discussed by [38,39]. They found that both annual and seasonal temperatures had significantly decreased at Gilgit, Gupis, Ushkore, and Yasin over the period 1991-2013.
A recent study [12] comprehensively discussed the Karakoram Anomaly in regard to balanced, or slightly positive, glacier budgets, resulting in advancing glacier termini and prevalent glacier surging activities. These conflicting signals are linked with negative trends in maximum and minimum temperature and increasing precipitation. Some recent studies have also endorsed the stability of glaciers in the Karakoram region, taking into account glacial mass balance, expanding SCA and hydroclimatic trends [6,13,35,38,65]. Our results are in partial agreement with those of abovementioned studies in verifying the increasing snow cover in the Gilgit basin, consistent with increased precipitation and significant decreases in temperature [38], however, snowmelt contribution dropped but glacier melt contribution increased during 2001-2012. Runoff within the GRB is mainly attributed to snowmelt, followed by glacier melt and rainfall [13]. Our findings, based on the results of SRM and SPHY, are consistent with the findings of [13], as we also noted the same contribution of discharge essentially depending upon the snowmelt, followed by the glacier melt and rainfall. We also noted the contribution of snowmelt runoff increased during 2001-2004 but dropped consistently after 2004 and continued up to 2012. On the other hand, we observed a consistently increasing trend of glacier melt contribution from 2004 onwards. According to some glaciological studies, eastern and central Karakoram glaciers in UIB are experiencing negative mass balances and consequently, are retreating [9]. Some recent studies in relation to glacier elevation changes [4,67,68] based on ICESat-1 and differential Global Positioning System (dGPS) data collection highlighted the role of debris cover in the Astore, Gilgit, and Hunza basins. They suggested a slight increase in the debris-covered area and less melting rate inversely correlating with the debris layer thickness [69]. The increased glacier contribution is inconsistent with less melting rate but likely to be attributable to monsoonal rains [37] (increased rainfall contribution) or associated with negative mass balances and rapid glacier retreat [9]. These findings still require energy-balance estimation and snow/ice albedo properties to establish a link with such complex anomaly. However, precisely assessed field-based DDF, debris-cover thickness and glacier surge activities may help to bridge the gap between model simulations (simplified parameterization) and in situ observations to explain such anomaly in a more refined way [12].
We found SRM is more sensitive to snow cover and temperature compared to precipitation. Specifically, the TLR is the most sensitive parameter; previous studies suggest −0.68 • C/100 m for the Karakoram Range. However, our updated analysis proposed −0.53 • C/100 m for the GRB. SRM being a conceptual model measures rainfall as snow (below T crit ) and vice versa. Therefore, decreased liquid precipitation at higher elevation zones (>5000) was observed throughout the study period (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). We observed that the GOF values for both models were efficient, but SPHY is more efficient in terms of NS and R 2 values, however, RMSE and bias were higher for both models. SRM's efficiency can be improved significantly by considering glacier melt contribution to runoff. The findings of the present study are also consistent with those of [29], as they suggested less model efficiency of SRM over the Hunza River Basin from July to August. We also support the notion of glacier and snowmelt contribution to river flows. SRM does not account for the glacier melt by default which causes error even when it is applied in the less glaciated region.
According to previous studies, a 1 • C rise in temperature will result in 33% and 16% increase in summer runoff of Hunza and Shyok rivers [14]. The river runoff response to climate is different for different basins largely depending upon the methods, hypothesis, and basin characteristics. Based on a linear regression analysis, [14] allocated the Karakoram region into high, middle, and low altitude based on glacier area, runoff dependence, and precipitation. According to these authors, the high altitude region (Shyok, Shigar, and Hunza) annual/summer runoff is attributed to mainly summer temperature, followed by middle (Astore, Gilgit, and Swat) and low altitude (Khan Khawer and Siran) in which flow is solely controlled by the winter precipitation (solid) and winter (liquid)/monsoon, respectively. This region can be further divided into the glacier dominant Karakoram Range (Shyok, Shigar, and Hunza) and snow dominant western Himalaya and Hindukush (Astore and Gilgit) [13]. The recent studies carried out in Gilgit and neighboring basins (Hunza and Astore), based on a similar hydrological modelling approach [28,29,37], suggested increasing flows between 1986 and 2010, attributed to increased temperature. Our updated analysis also suggested increased future flows attributed to increasing temperature. We observed that the flows were not increasing according to the temperature increase as there would be less snowmelt and glacier melt due to intensified warming, so overall water availability in terms of cryospheric extent will likely decrease during 2070-2100. These results are consistent with [37], who reported robust decreases in future runoff during the same period within the Hindukush mountainous range and Himalayan region [70] during the 21st century. The increases in temperature trends are in the phase with global warming, resulting in the negative mass balance of Himalayan and Hindukush glaciers [9].
The present study inclines more towards the uncertainty caused due to the selection of modelling parameters-DDF particularly-in terms of debris-covered and clean-ice glaciers rather than the propagation of uncertainty caused by future climate and hydrology. Some earlier studies emphasized that future climate in the upper Indus basin is highly uncertain as none of the current state-of-the-art GCMs and RCMs satisfactory simulate the monsoon and westerly dynamics in the region, making the reliability of future scenarios questionable [71][72][73]. The Gilgit sub-basin is subject to large uncertainties in both RCPs (4.5 and 8.5) resulting from uncertainty in the predictability of the monsoon, autumn precipitation and summer air temperatures [37]. Recognizing this, we employed an ensemble of a limited number of GCMs selected on the basis of previous research. While the models represent only a sample of the total uncertainty space, the results were quite different in terms of magnitude during various data series depending upon the scenarios. We found significant increasing temperature for both GCMs during the second data series.
If realized, increased summer flows (peak-water) in the region could prove beneficial for a range of sectors like hydropower and irrigation, but only over the short to medium term and if not associated with extreme events. Importantly for a water management perspective, the long-term projections indicate declining water resources in the region in terms of snow and glacier melt. We concluded that over the long term, that glacier and snow meltwater will be insufficient to need the needs of increasing demands associated with increased human activities and/or population. Water management facilities in Pakistan currently lack the capacity to store sufficient water resources during the peak-water season (ablation) for use during drought/low flow periods. The sensitivity of the water sector was recently highlighted when the storage capacity of the Tarbella, Mangla, and Chashma hydropower reservoirs was compromised, due to heavy silt loads (Tarbella and Mangla), resulting in dam operations running on much reduced capacity. To cope with limitations in storage capacity, increasing existing, or the development of new, storage capacity, as employed at the Mangla reservoir, will be required. Moreover, being an agricultural-based country, crop irrigation is heavily reliant on the water storage capacity within the Indus basin. The inability to store water is particularly highlighted by the irrigation requirement of the Rabi crop, which grows from October to March. In order to prevent structural damage during this period, dams must release excess water, which results in severe water shortage during the Rabi growth season. If climate change increases the seasonality of flows, existing infrastructure will not be capable of storing sufficient water resources for use during the low flow periods and result in increased vulnerability of the region to climate change.

Conclusions
Our study suggests that SRM simulates flows reasonably well during April and May, but during the peak ablation period, its performance diminishes when glacier melt contribution typically dominates snowmelt runoff, particularly evident during July and August. SRM is more sensitive to temperature and DDF than precipitation at high altitudes. The present study finds a marginally increasing, more stable, snowpack during the period 2001-2018, consistent with enhanced winter precipitation and decreased summer temperature within the Gilgit basin. Increasing precipitation is likely associated with increased snow cover at high altitude across the region. Our results revealed that the most sensitive parameter in the simulation of runoff is the lapse rate and DDF for debris-covered and debris-free glaciers, which dictate the shape of the hydrographs of observed and measured runoff based on over and underestimation of runoff during the ablation season. Model efficiency can be improved further by considering the climatic data (temperature and precipitation) obtained at high altitude weather stations. Integrated mountain water resource management could be carried out more efficiently by applying this model in other snow-fed basins considering field data collected through physical glacier melt observations, which would allow for greater intercomparison of satellite and field-based data over UIB. Our analysis indicates that SRM was capable of simulating daily flows during the summer months and may provide a suitable tool for use in water resources engineering and management.
In general, the efficiency of the hydrological models employed could be increased using field-based data collected through frequent visits, particularly during early and late ablation periods. The results should also be validated through glacier mass balance analysis based on physical data collection from accessible glaciers. With regards to SPHY, efficiency could be further improved using field-based melting rates of debris-covered and clean-ice glaciers during early and late ablation periods. Moreover, high altitude data collection of precipitation and temperature is as important as glacier mass balance data and should be prioritized as part of any future monitoring program. With improved monitoring and subsequent model efficiency, more robust findings would further underpin the design of improved sustainable water management policy under changing climate.
The modelling approaches employed in this study also highlight the need for research to focus on improved mapping, particularly in terms of permafrost extent and ice volumes. These two entities still need high-resolution classification/mapping based on the individual basin-scale hydrology for more precise runoff assessment in a snow/glacier-fed basin. While new satellite technology will assist in this, the need for ground-based measurements remains. The Karakoram anomaly might have positive impacts, in terms of positive basin storage; however, its persistence into the future remains uncertain. It is essential that long term monitoring programs be developed, or existing programs enhanced and expanded, to ensure that scientists and water managers have sufficient information to develop new or adapt existing water management plans in a timely manner.