Modeling Hydrological Response to Climate Change in a Data-Scarce Glacierized High Mountain Astore Basin Using a Fully Distributed TOPKAPI Model

Abstract

Water scarcity is influencing environmental and socio-economic development on a global scale. Pakistan is ranked third among the countries facing water scarcity. This situation is currently generating intra-provincial water disputes and could lead to transboundary water conflicts. This study assessed the future water resources of Astore basin under representative concentration pathways (RCP) 4.5 and 8.5 scenarios using fully distributed TOPographic Kinematic APproximation and Integration (TOPKAPI) model. TOPKAPI model was calibrated and validated over five years from 1999–2003 with a Nash coefficient ranging from 0.93–0.97. Towards the end of the 21^st century, the air temperature of Astore will increase by 3°C and 9.6 °C under the RCP4.5 and 8.5 scenarios, respectively. The rise in air temperature can decrease the snow cover with Mann Kendall trend of −0.12%/yr and −0.39%/yr (p ≥ 0.05) while annual discharge projected to be increased 11% (p ≤ 0.05) and 37% (p ≥ 0.05) under RCP4.5 and RCP8.5, respectively. Moreover, the Astore basin showed a different pattern of seasonal shifts, as surface runoff in summer monsoon season declined further due to a reduction in precipitation. In the spring season, the earlier onset of snow and glacier melting increased the runoff due to high temperature, regardless of the decreasing trend of precipitation. This increased surface runoff from snow/glacier melt of Upper Indus Basin (UIB) can potentially be utilized to develop water policy and planning new water harvesting and storage structures, to reduce the risk of flooding.

1. Introduction

The climate change impact on the glaciated Hindu Kush–Karakoram–Himalaya (HKH) mountain region is the subject of extensive debate for many years and still unresolved due to the scarcity of ground station data especially for high altitude [1]. This is further complicated due to poor understanding of the interaction between summer South-Asian monsoon and winter westerly disturbance systems, which leads to the disagreement between global, regional, and local-scale climate change projection [2]. Many studies have reported extensive changes in the spatio–temporal distributions of precipitations (rainfall and snowfall) [3]. Since air temperature determined whether the precipitation will fall as a solid (snow) or liquid (rain), therefore the timing of the precipitation received by a basin will determine its amount, type, and ultimately the availability of water downstream [4].

The Indus is the backbone for the economy of Pakistan, particularly for Punjab and Sindh because these provinces are responsible for most of the nation’s agricultural production. Also, the water stored in its reservoirs is used to generate electricity and provides the main supply of drinking water in Pakistan [5]. More than 90% of the Indus River lowland flow is potentially contributed by snow/ glacial melt from the UIB whereas 60% of the total annual runoff to the Indus River is attributed by the melting of perpetual snow and glaciers which exist above 3500 m in the sub-basins of UIB including Shyok, Shigar, Astore, Gilgit, and Hunza. Since the economy of Pakistan and a large number of the population around 207 million people [6] are highly dependent on snow and glacial meltwater of UIB, so it is essential to assess the response of these water resources under IPCC climate change scenarios.

Many studies [7,8,9,10,11] have analyzed the impact of climate change on water resources of the UIB using hydrological models. However, most of these hydrological models are simple lumped or conceptual, and simple degree-day models, i.e., Snowmelt Runoff Model (SRM), Hydrologiska Byråns Vattenbalansavdelning (HBV), Soil and Water Assessment Tool (SWAT), and water and energy. budget-based distributed hydrological model WEB-DHM-S model rather than fully physical [12,13,14] and spatially distributed [14,15]. Fully physical spatially distributed hydrological models solve the numerical equations related to conservation of momentum, mass, and energy both in time and space [16] and hence, require an extensive input of several datasets. Whereas conceptual hydrological models only approximate the general hydrologic processes of an area [17], and maybe less challenging regarding data inputs. The contrasting hydro-meteorological regimes of UIB are highly sensitive to an abrupt and changing climate, mainly due to changing precipitation (rainfall and snow) and temperature variability [18,19]. Due to the climatic variability, scarcity of long-term and consistent data, the resolution and accuracy of Global Circulation Model outputs and less understanding about the response of glaciers and snow cover changes in the UIB, it is very challenging to model the accurate hydrological response of these catchments to climate change [20,21,22,23]. Archer [24] indicates that a 1 °C increase in average summer temperature would result in a 16% and 17% increase in summer runoff due to glacial melt in the Hunza and Shyok sub-basins, respectively. Another study conducted in the same basins by Fowler and Archer [25] through trend analysis of meteorological datasets reported that 1 °C rise in average annual temperature would accelerate glacial melt. Singh and Bengtsson [26] applied a lumped snowmelt model (SNOWMOD) on Sutlej river basin and concluded that an increase of 1–3 °C temperature would accelerate the glacier melt by 16%–50%. Usually, the researchers did not relate the glacier dynamics with the hydrological cycle like evapotranspiration, surface runoff, and snow water equivalent through hydrological modeling [27]. Most of the studies in UIB applied simple temperature index-based methods [28] and considered possible depletion in glacier dynamics [12,19,26,28]. Though these research studies give some crucial understandings about the impacts of changing climate on water resources of UIB but, they did not incorporate the possible future changes in glaciers, since 60–70 percent of the flow is generated in the UIB due to glaciers and snowmelt as well as spatial dynamics of other hydrometeorological variables at sub-basin level.

Based on above knowledge gaps present in the hydrological studies of UIB there is dire need to assess the impact of changing climate on the hydro-glacial sub-basins with a model which incorporates the physical process and spatial variability of both glaciers and hydrometeorological variables. The most appropriate approach to resolve this issue is the implication of a fully-distributed physically-based hydrological model that can simulate the physics of the main governing processes at the sub-basin level of UIB [29]. Since Freeze and Harlan [30] introduced the first physically based spatially distributed hydrological model, since then many models have been developed like MIKE Système Hydrologique Européen (SHE) [31], the Institute of Hydrology Distributed (IHD) model [32], Coupled Routing and Excess STorage (CREST) [33], TOPKAPI (TOPographic Kinematic APproximation and Integration) [34,35] and many others. The TOPKAPI model could be used in a data-scarce catchment because it does not need calibration of the parameters while only incorporation of watershed characteristics data will produce highly reliable results as compared to other models like conceptually based model, which need observed hydrograph for automatic calibration which reduces its application in ungauged basins [36]. Another advantage of this model is that it incorporates the glacier boundaries to compute the point scale glacial discharge. There are several other models that incorporate the glacier boundary like Soil and Water Assessment Tool (SWAT), Spatial Processes in HYdrology SPHY and Hydrologiska Byra ˚ns Vattenbalansavdelning (HBV) model but they simulate the glacier melt based on simple Temperature Index models (TIMs) or Degree-Day Models (DDMs) which relies on empirical relationships between air temperatures and melt rate [28,37]. TOPKAPI glacial and snowmelt component is based on the Enhanced Temperature Index (ETI) model [38], which incorporates the information of shortwave radiation and albedo along with air temperature to give the more reliable results as compared to TIMs and DDMs. In previous studies, the applications of such hydrological models in UIB were very limited due to their requirements of extensive spatial data inputs [1,39]. The UIB region is well known for the scarcity of hydrometeorological and particularly glacier mass balance data [40,41]. One or two meteorological stations and one hydrological station is present in each sub-basin of UIB which gives limited and sparse spatial coverage of in situ observations in the area. Moreover, there is no high altitude (>4500 m) hydro-meteorological station located in the sub-basin of the UIB. Recently, remotely sensed satellite data have emerged as a viable option to provide datasets to land/hydrologic models or addition to in situ observations because of their spatiotemporal availability over substantial ungauged regions like UIB. Several hydro-meteorological parameters can be derived from remote sensing-based satellite data [42,43,44,45,46,47,48]. In this study, we have also incorporated satellite remote sensing data in the TOPKAPI model to assess the impact of future climate changes on the hydrology of the glaciated Astore river basin. These previous studies did not incorporate the spatio–temporal variability of the basin glaciers in a fully distributed hydrological model to better understand the hydro-meteorological behavior under the looming future climate change scenarios. The specific objectives of the study were 1) calibration and validation of fully physical distributed TOPKAPI hydrological model to simulate the surface runoff in the snow and glacierized Astore river basin daily; 2) to project the impact of future climate change on the hydrology of Astore basin.

2. Materials and Methods

2.1. Study Area

The Astore catchment (a sub-basin of the UIB) located at 35° 33′ N, 74° 42′ E, was selected (Figure 1) for hydro-glaciological modeling under IPCC climate change scenarios. The Astore River basin is located in a northwestern Himalayan arc, with an elevation range 1198 m to 8069 m and total drainage area is about 3988 km², Figure 2 shows the hypsometric curve. About ~543 km² area is covered with glaciers with the highest of the peak of Nanga Parbat, which is ranked ninth highest mountain in the world [49]. The average annual temperature during 1998–2012 was 9.9 °C at the station located in valley Rattu (2718 m) and −2.9 °C at the higher altitude Rama station (3179 m).

2.2. TOPKAPI Model

TOPKAPI is a fully distributed physically-based hydrological model that can simulate the high spatial and temporal information on the hydrological state of both small and large watershed [14] and mountainous catchments [50]. The primary inputs required by TOPKAPI model are terrain data or Digital Elevation Models (DEM), soil maps, land cover, glacier boundaries, meteorological data, i.e., temperature and precipitation. Several inputs of soil and hydrological control parameters are required to run the model. Soil parameters includes L (soil depth in m), k_s (vertical and horizontal permeability at saturation), ϑ_r (residual water content), ϑ_s (saturated water content), α_s (non-linear soil exponent of the transmissivity) whereas, routing parameters includes n_o (manning coefficient for the land cover in m^−1/3 s⁻¹) and n_c (manning coefficient for the channel in m^{−1 /3} s⁻¹). Fully distributed physically-based models provide high resolution and more reliable simulations of hydrological behavior towards climate change than empirical or lumped models, which heavily depends on the parameter’s calibration. These models require extensive data input and need to be calibrated and validated to generate all relevant processes accurately [50].

TOPKAPI model has been modified extensively for various application in the mountainous regions to simulate water balance at catchment and sub-catchment scales. The model generates high spatial and temporal outputs for every cell of entire catchment which include time series of hydrologic quantities, i.e., precipitation, air temperature, potential and actual evapotranspiration, snow cover, percolation, surface flow, channel flow, soil moisture and glacier /snow water equivalent. The TOPKAPI model uses the Enhanced Temperature Index (ETI) approach to simulate the glacier and snowmelt component. ETI is an intermediate complexity step between an empirical TIMs and an Energy Balance Models (EBMs) [28,38]. ETI model computes the melt component when the air temperature at specific location of grid cell iT_i (°C) exceed the threshold air temperature T_T.

M_i= TF.T_i + SRF.I_Gi * (1 - α_i)

(1)

where M_i is the melting in a particular cell i (mm water equivalent (w.e.) per day d⁻¹, TF is the temperature factor (mm w.e. K⁻¹ d⁻¹), SRF is the shortwave radiation factor (mm w.e. m² W⁻¹ d⁻¹), I_Gi is the incoming shortwave radiation in (Wm⁻²), and α_i is the surface albedo.

The advantage of ETI is that data requirement is not much higher as compared to EBMs. In addition to the air temperature term, the ETI model includes a shortwave radiation term that incorporates incoming solar radiation and albedo that help to account the internal energy of the glacier to provide reliable results. So the TOPKAPI model proved to be more efficient to simulate the glacier melt than the other TIMs and DDMs. With several advantages of TOPKAPI model, it has some limitations that it does not provide the separate output of changes in glacier dynamics and gives the combined hydrograph from glacier melt, snowmelt and precipitation. These output data can be useful to analyze the behavior of the entire basin/sub-basin during severe hydro-meteorological events and can be utilized for preparedness, mitigation, and relief in case of disaster.

2.3. Data Input to Model

The TOPKAPI model input data includes the Digital Elevation Model (DEM) with a spatial resolution of 30 × 30 meters grid, the soil data, the land cover raster, and the hydro-meteorological data (daily average temperature, precipitation, and streamflow).

Figure 3 shows the hydrological and topographical variables derived from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global DEM in the ArcHydro module of Esri’s ArcGIS [43].

Landsat 7 Enhanced Thematic Mapper Plus (ETM+) with a spatial resolution of 30 m was used to generate a Landover map of the Astore basin. Preprocessing (i.e., geometric and atmospheric correction) of Landsat image was done using different filters followed by application of the Maximum Likelihood Classifier (MLC) technique to classify the land cover features that depend on the second-order measurements of a Gaussian probability density function model [51,52]. Figure 4 shows the major land cover classes of Astore basin. TOPKAPI model requires the manning roughness coefficient and crop factor value of all months for each land-use class.

Glacier boundaries were delineated using a morphometric glacier mapping (MGM) technique. The basis in MGM is the use of topographic thresholds derived from a DEM and the multispectral data along with ratio images for the automatic delineation of glacier boundaries. The filters were generated by threshold values defined for the extraction of glaciers. The Normalized Difference Snow Index (NDSI), was used to delineate the glaciers (Figure 4). To further refine the glacier’s boundaries, other filters like Normalized Difference Vegetation Index (NDVI), Land Water Mask (LWM), mean hue [46], slope, and elevation threshold values were used. The delineated glacier outlines were validated using Google Earth’s high-resolution imagery [53].

The FAO global soil dataset [47] with a spatial resolution of 30 arc-second was used to generate the soil data of the Astore basin (Figure 5). The TOPKAPI model requires soil property data such as the soil depth, horizontal and vertical soil hydraulic conductivity, available moisture content, and texture. These soil hydraulic parameters were derived using the PedoTransfer Function (PTF)-based software SoilPAR 2.0 [48] and Soil Water Infiltration and Movement (SWIM) model [49].

Meteorological data of Rama and Rattu stations (1999–2003) comprising of daily average temperature and precipitation and daily flow records (1980 and 2010) at Doiyan gauging station (1583 m asl) of the Astore River was input in the model. This hydro-meteorological dataset was used for model calibration (1999–2000) and validation of model (2001–2003).

Future Climate Projection Data

For future climate change modeling, the Asian Precipitation Highly Resolved Observational Data Integration Toward Evaluation of Water Resources (APHRODITE) dataset was used that has been validated through the actual meteorological stations. The APHRODITE data products V1204R1 and V1101R2 were downloaded for daily temperature and precipitation, respectively. The APHRODITE data was selected to downscale the climate data that can provide the high spatial resolution of 0.1° and 24 h temporal output to fulfill the data input requirement of fully distributed TOPKAPI model. However, APHRODITE data have the limitation that the precipitation measurements are underestimated at a high altitude of UIB, so the bias-corrected APHRODITE temperature and precipitation data were incorporated in the model at a spatial resolution of 0.1°.

In addition to historical data of APHRODITE, Global Climate Models (GCMs) projections generated from the Coupled Model Intercomparison Project Phase 5 (CMIP5) under Representative Concentration Pathways (RCPs) scenarios was also incorporated to assess the climate change impacts on the future hydrology of Astore basin. Four GCMs [54] (i.e., CCSM4, CanESM2, HadGEM2–ES, and GFDL–ESM2M) were selected based on the following criteria (a) models that showed strong correlation (≥ to 0.8) with baseline APHRODITE data (1960–2010), (b) Normalized Root Mean Squared Error (NRMSE) ≤ 0.1 and (c) Normalized Standard Deviations (NSD) was within the range of ±0.4 to that of NSD of APHRODITE time series data.

CCSM4 data with very high correlation (0.94) and less NRMSE value (0.01) was biased corrected under the method of Linear Interpolation and Bias Correction (LIBC) technique [55,56,57]. Linear Interpolation and Bias Correction Method (LIBC) utilizes station data and high-resolution topography, and it is the combination of the geospatial Linear Interpolation (LI) technique in relation with Bias Correction and Spatial Disaggregation (BCSD) technique. BCSD takes account of a correction factor to generate the biased corrected output of GCMs. Bias corrected climate parameters (i.e., temperature and precipitation) of APHRODITE were chosen and interpolated at a resolution of 0.1 degrees for the baseline period 1975–2010. The daily dataset was converted into monthly time series and manipulated to get the climatology and standard deviation of the parameter. A similar process was done for GCMs data, first linearly interpolated at the resolution of APHRODITE, and monthly time series data was manipulated to get the climatology and standard deviation of the parameter. After that correction factors were calculated by the ratio of climatological values of observations to the reference GCMs. Bias corrected APHRODITE datasets have been used in downscaling of baseline GCMs to generate future scenarios daily from the year 2020–2100 with the spatial resolution of 0.1°. This technique eliminates various inaccuracies and the biases related to GCMs data. The result was high resolution gridded projections of temperature and precipitation that used effectively in the TOPKAPI model to predict the hydrological response of the study area at a spatial resolution of 15m grid daily.

2.4. Sensitivity Analysis of the Model Parameters

The sensitivity analysis of various parameters of a hydrological model is assessed to identify which parameter has the highest impact on the model output. The sensitivity of the TOPKAPI model parameters was determined by changing their default values to see each of the parameter’s sensitivity. The initial default values of selected parameters were varied within the upper and lower limits as suggested by literature and documentation of the TOPKAPI model [58,59] as highlighted in

Table 1. Values of the TOPKAPI (TOPographic Kinematic APproximation and Integration) model parameters estimated a priori from data and literature, and values obtained from the calibration procedure.

Parameter	Pre-Calibration Value	Origin and References	Post-Calibration Value
Cell specific
Ground Slope tan, β	1.7E-4–1.81E-1	[62]
Channel Slope tan, βc	1.7E-4–1.81E-1	[62]
Depth of surface soil layer (m), L	0.1–0.81	[67]	0.1–4.5
Saturated hydraulic conductivity (m.s−1) Ks	1.09E-8–1.09E-6	[72]	1E-12–1E-002
Residual soil moisture content (cm3/cm3) θr	0.01–0.047	[14]
Saturated soil moisture content (cm3/cm3) θs	0.1–0.46	[14]
Manning’s surface roughness coeff. (m−1/3.s−1) no	0.03–0.1	[61,73]	0.001–0.25
Manning’s channel roughness coeff. (m−1/3s−1) nc	0.02–0.07	[15]	0.11–0.43
Non-linear soil exponent αs	2.5	[15]
Crop factor kc	1.0	[62]	0.01–1
Horizontal dimension of cell (m) X	1 000	[74]
Max. channel width at outlet (m) Wmax	10	[62]
Min. channel width (m) Wmin	1	-
Area required to initiate channel (m2) Athreshold	25,000,000	[51]

. The model’s sensitivity against these parameters was evaluated by comparing the observed and simulated discharge through Nash efficiency statistics. Based on the visual and statistical agreement between the observed and simulated discharge the most sensitive parameters, affecting the hydrology of the basin were selected [60] for tuning in the calibration process of the TOPKAPI model. In this study, the most sensitive parameters were the depth of surface soil layer (L), saturated hydraulic conductivity (K_s), Residual soil moisture content (θ_r) and saturated soil moisture content (θ_r).

Model Calibration and Validation

The TOPKAPI model, like other physically-distributed models, is subjected to several uncertainties associated with the input data. These uncertainties may be due to topographic variation, land cover types, soil properties, approximation methods and tables used to infer physical parameters from the data, and the approximations introduced by the scale of the parameter representations [61]. For these reasons, Liu and Todini [15] suggested that the calibration of parameters is still necessary, but all studies dealing with the TOPKAPI model suggested that the model calibration can be achieved through simple hit and trial method, as we had also done in this current study for Astore basin.

In this current study, an automatic model calibration method i.e., Ordered Physics-based Parameter Adjustment Method (OPPA), as proposed by Vieux, et al. [62] was used for the model calibration. The OPPA is the most effective and time-saving method recommended for the calibration of physically distributed hydrological models, as this method calibrates a model in a specific order [62,63,64]. Initially, the method adjusts the runoff generation parameters to minimize the discharge volume objective function and then those parameters which controlled the runoff routing to minimize the discharge timing objective function. Through the sensitivity analysis and literature Liu, Martina and Todini [14] it was found that soil depth L and soil conductivity K_s are the most sensitive runoff generation factors while the runoff routing is mainly controlled by the Manning roughness of the channel n_c and the over-land surface.

Another important parameter, i.e., soil moisture V_{s_initial} was also tuned because of its high sensitivity towards the runoff generation. To initialize the soil moisture parameter in the model, it was assumed that there was no surface water over the slopes and the water depth in a generic channel cell increases as a linear function of the channel width of the cell. At the beginning of the calibration process, the initial value for soil moisture percentage was set according to dry condition. Then the model was run to find the soil moisture situation at a specific time when the soil should be in the same state as the initial condition. This OPPA method provides accurate and realistic values based on the values for the initial condition, especially concerning the spatial distribution of the saturation percentage. According to Vieux, Cui, and Gaur [62] and others [65,66,67], the hydrological models cannot be calibrated at each pixel level. However, a multiplicative factor can be used for the whole basin. In this study, the four multiplicative factors L (for the soil depth), K_s (for the hydraulic conductivity), n_o (for the overland roughness) and n_c (for the channel roughness) were applied for the model calibration.

The parameters were grouped into (L, K_s, V_s) and (n_o, n_c) because they were independent of each other, i.e., the change in the first group was affecting the quantity and change in second was affecting the timing of runoff at Astore basin (Table 1). Three different statistical measures as recommended through literature, i.e., Index of Agreement (d), R² and Nash–Sutcliffe Efficiency (NSE) coefficients were used to assess the efficiency of calibration by comparing modeled and observed discharge volumes aggregated on a daily and monthly basis [68,69,70,71].

Before further application of a hydrological model for future forecasting, it is recommended that the validation of the model should be done for a better understanding of spatio–temporal performance of the model. Hence, after the model calibration, validation of TOPKAPI model was done on the remaining hydrological data (2001–2002) to evaluate the performance of simulation and prediction model based on the statistical value of R², Index of agreement between the simulated and observed discharge and Nash–Sutcliffe Efficiency (NSE) coefficients [75].

3. Results and Discussion

3.1. Model Calibration

The TOPKAPI model calibration (1999–2000) had a significant (p ≤ 0.05) relationship between observed and simulated Astore River discharge with an R² value of 0.96 and Nash–Sutcliffe efficiency value of 0.95 (

Table 2. Statistical evaluation of modeled simulations during the calibration period 1999–2000.

Event	R2	NSE	d1
1999	0.97 *	0.95 *	0.93 *
2000	0.93 *	0.94 *	0.97 *
1999 to 2001	0.96 *	0.95 *	0.95 *

* indicate significance at α = 0.05 and d is the Index of Agreement [76].

). Figure 6 shows the comparison between the observed and TOPKAPI model simulated discharge along with temperature and precipitation in the basin.

The maximum average flow during the calibration period was observed throughout the summer season (model 188m³.s⁻¹, observed 190 m³.s⁻¹), whereas discharge started increasing in late April and reached its peak in June. This hydrological response of the Astore basin could be due to summer monsoonal liquid precipitation in the basin. Additionally, the higher temperature 12 °C than the average temperature in the basin during the summer season caused increased melting of snow/glaciers resulting in peak discharges of Astore basin. Several research studies also reported a similar trend in hydro-meteorological data [77,78,79]. Whereas, moderate river flow was observed in the late summer (model 81 m³.s⁻¹, observed 85 m³.s⁻¹) with a decrease in discharge started in late August till mid of October. This hydrological pattern of the basin is also in confirmation with temperature trends during a similar period. In winter, the river discharge was very low (model 30 m³.s⁻¹, observed 38 m³.s⁻¹), which shows a close relationship with the temperature of the basin. As the temperature of the basin is at its minimum below 0°C, hence the precipitation in winter is solid, i.e., snowfall. This snowfall caused the accumulation of glaciers and snow in the high altitudes (>4500 m) [52,79,80].

3.2. Model Validation

Figure 7 shows a comparison of simulated and observed discharge at Doiyan (Astore River Basin outlet) along with the trend of the temperature and precipitation, during the validation period 2001–2003.

Table 3. Statistical summary of flow simulation during the validation period 2001–2003.

Event	R2	NSE	D
2001 2002 2003	0.98 *	0.95 *	0.94 *
	0.97 *	0.92 *	0.94 *
	0.97 *	0.93 *	0.96 *
2001 to 2003	0.96 *	0.93 *	0.95 *

* indicate significance at α = 0.05 and d is the Index of Agreement [76].

depicts the R², Nash–Sutcliffe efficiency (NSE), and index of agreement (d) statistics. The 2001 validation period had higher R² and NSE values. However, there was no significant difference.

The maximum average flow during the validation period was observed throughout the summer season (model 199m³.s⁻¹, observed 257m³.s⁻¹, whereas discharge started increasing in late April and reached its peak in June during 2001–2003. This hydrological response of the Astore basin could be because of the monsoonal liquid precipitation in the basin. Additionally, the increase in temperature during the season caused an increased melting of snow/glaciers resulting in the peak discharges. Whereas, moderate river flow was observed in the late summer (model 98 m³.s⁻¹, observed 84 m³.s⁻¹) with a decrease in discharge started in late August till mid of October. This hydrological pattern of the basin is also in resemblance with temperature trend for the same months. In winter, the river discharge was very low (model 28 m³.s⁻¹, observed 39 m³.s⁻¹), which shows a close relationship with the temperature of the basin.

It was observed that about 89.74% of runoff volume generated in the summer season (May–September) due to higher glacier/snow meltwater and summer monsoonal rainfall, whereas insignificant runoff, was observed during winter since most precipitation received as snow in the basin during 1999–2003. Winter runoff was mainly base-flow generated from snowmelt water that was stored in the soil.

Figure 8 shows significant relationships between observed and simulated Astore river daily discharge for the calibration and validation period. The minimum NSE coefficient value never dropped below 0.93 during the entire validation period. The best performance of the model was found over the calibration period with the NSE coefficient value of 0.95, and according to Moriasi, et al. [81], a hydrological model is well calibrated for a runoff if monthly NSE ≥ 0.65. Besides the runoff generation, the TOPKAPI model also generates several other outputs as given in

Table 4. Average monthly flow statistics of flow simulation during the validation period 2001–2003.

Month	Temperature (ºC)	Precipitation (mm)	Soil Moisture (%)	Snow Cover (%)	SWE (mm)	ETP (mm)	ETA (mm)	Percolation (mm)
Jan	−11.33	1.57	17.02	91.35	229.84	0.00	0.00	0.00
Feb	−10.31	3.09	16.09	91.76	290.47	0.00	0.00	0.00
Mar	−6.77	4.56	15.82	89.33	404.26	0.00	0.01	0.00
Apr	−1.00	4.70	18.55	89.76	453.86	0.00	0.11	0.04
May	4.29	2.27	27.46	61.41	315.41	0.52	0.33	0.35
Jun	8.16	1.09	31.18	31.00	88.21	1.08	0.46	0.27
Jul	11.00	0.86	27.63	10.62	18.16	1.46	0.47	0.03
Aug	10.71	1.00	25.05	8.17	12.02	1.39	0.44	0.02
Sep	6.39	1.32	22.99	11.40	11.98	0.81	0.31	0.02
Oct	0.17	1.09	20.45	41.39	18.63	0.06	0.08	0.02
Nov	−4.95	4.09	19.00	72.19	71.83	0.00	0.01	0.01
Dec	−9.46	3.13	17.82	86.99	187.91	0.00	0.00	0.00

.

An analysis of snow cover, percolation, evapotranspiration, soil moisture, temperature, snow water equivalent (SWE) and precipitation are vital for understanding the integrated hydrologic response of a catchment. The accurate estimation of basin-wide SWE is the most important unsolved problem in mountain hydrology [82], particularly in the HKH region. The SWE is the amount of water contained within the glacier, which can be taken as the amount of water that would theoretically result if the entire snowpack is being melted instantaneously [83]. TOPKAPI model efficiently simulated the SWE in the Astore river basin, and the results showed that the maximum SWE was found in April (453.86 mm) just after the end of the snow accumulation period, i.e., winter and the lowest value was found in September (11.98 mm) at the end of the summer period. The other significant output of the TOPKAPI model is snow cover extent; the results indicate that the average snow-covered in two seasons, i.e., in winter (November to April) and summer (May to October) was 91 % and 11% respectively. A similar snow cover trend was reported by Atif, Iqbal, and Mahboob [78], they reported that the snow cover in the Astore basin varies from 92% during winter to 15% in summer seasons.

3.3. Future Climate Projections (2010–2100)

Very fine resolution future climatic data (temperature and precipitation) acquired by the CMIP5 model were downscaled under the two Representative Concentration Pathways scenarios (RCP4.5 and RCP8.5) and used to forecast the future hydrology of Astore. Monthly climate change for the three periods (2010–2040, 2041–2070, and 2071– 2100) is summarized in

Table 5. Monthly average temperature variation (°C) under near (2010–2040), middle (2041–2070), and future (2070–2100) climate change scenarios.

Period	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	Annual
RCP4.5
2010–2040	−10.4	−8	−3.0	3.2	8.8	12.1	15	14.9	10.8	4.9	−1.2	−7.3	3.4
2041–2070	−6.3	−7.8	−5.4	0.4	7.0	11.9	14.8	16.7	14.9	10.5	6.2	0.0	5.3
2070–2100	−2.1	−5.7	−5.9	−3.4	2.7	8.5	12.2	15.6	17.4	14.6	9.8	3.9	5.7
2010–2100	−6.2	−7.1	−4.8	0.0	6.1	10.8	14	15.7	14.4	10	5.0	−1.1	4.8
RCP8.5
2010–2040	−10	−8.4	−2.8	4.0	9.1	12.5	15.4	14.7	10.6	5.1	−1.3	−7.3	3.5
2041–2070	−3.6	−6.1	−4.2	2.3	8.4	12.2	15.4	18.3	17.0	12.6	8.1	2.2	6.9
2070–2100	5.3	1.0	−0.6	1.7	8.4	14.4	17.4	19.9	22.3	19.8	14.4	9.6	11.2
2010–2100	−2.8	−4.5	−2.6	2.7	8.6	13	16.0	17.6	16.6	12.5	7.1	1.5	7.2

and

Table 6. Monthly average precipitation variation (mm) under near (2010–2040), middle (2041–2070) and future (2070–2100) climate change scenarios.

Period	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	Annual
RCP4.5
2010–2040	1.3	1.8	2.0	2.0	2.0	1.9	1.6	1.7	1.4	0.9	0.9	1.1	1.6
2041–2070	1.3	1.5	1.4	1.6	1.9	1.6	1.8	2.2	1.8	1.0	1.0	1.2	1.5
2070–2100	0.9	1.9	2.6	2.2	2.2	2.6	2	1.5	1.6	1.1	0.9	0.7	1.7
2010–2100	1.3	1.8	2.0	2.0	2.0	1.9	1.6	1.7	1.4	0.9	0.9	1.1	1.6
RCP8.5
2010–2040	1.8	2.1	2.1	2.6	2.3	1.3	1.2	1.5	0.9	0.8	1.4	1.4	1.6
2041–2070	2.7	2.5	1.9	1.7	2.6	2.4	1.5	1.5	1	0.6	0.9	1.8	1.8
2070–2100	1.3	2.3	2.3	1.6	1.2	1.5	1.5	1.5	1.8	0.9	0.9	1.1	1.5
2010–2100	1.9	2.3	2.1	2.0	2.0	1.8	1.4	1.5	1.2	0.8	1.1	1.4	1.6

.

Towards the end of this century, an increase of 0.03 °C and 0.12 °C in annual temperature was observed under RCP4.5 and RCP8.5, respectively. The analysis showed that under RCPs the summer temperatures would be high as compared to the base data, i.e., 3.4 °C high in RCP4.5 and 6.4 °C in RCP8.5 (Figure 9). Research studies conducted by Kazmi has also supported our results. Rajbhandari, et al. [84] and Ul Islam, et al. [85], which concluded an increase in future temperature by 4.8 °C at the end of this century using PRECIS Regional Climate Model (RCM) model.

Under the RCP8.5 and RCP4.5 mean precipitation projection shows an increasing trend. Generally, the RCP8.5 scenario shows an 11% significant (p ≤ 0.05) increase at the end of this century. Several studies reported a similar increasing trend in precipitation in the UIB, for example, Beven [13] reported a 25% precipitation increase compared to the baseline period (2046–2065) based on the results of five General Circulation Models (GCMs). Also, Rajbhandari, Shrestha, Kulkarni, Patwardhan, and Bajracharya [84] forecasted an increase in precipitation up to 21% using the PRECIS fine-resolution RCM. Another study conducted by Forsythe, et al. [86] also reported an increase of precipitation in the UIB by applying a stochastic rainfall model. All of the above studies projected an increasing trend in the UIB precipitation regardless of using different methodologies to conclude. However, uncertainties exist in the magnitude of precipitation projected by these studies. Since we also found similar uncertainties regarding our precipitation projections, therefore it is important to conduct further studies at the basin level to reduce those uncertainties.

3.4. Future (2010–2100) Hydrological Flows Trend

The TOPKAPI model simulated a negative snow cover projections trend of −0.12%.yr⁻¹ and −0.39%.yr⁻¹ under IPCC climate scenarios RCPs 4.5 and 8.5, respectively during 2010–2100 (Figure 10). The corresponding decrease in snow cover ranged from 47% to 35% and 55% to 21% for RCP4.5 and RCP8.5 scenarios, respectively during the end of this century 2100. The mean snow cover area was 46% during the near future (2010–2040); 38% during the middle future (2041–2070) and 22% during far future projection (2071–2100) under RCP8.5. Whereas under RCP4.5, mean snow cover area was found to be 47% during (2010–2040); 39% during the middle future (2041–2070) and snow cover extent remained the same during (2071–2100).

Figure 11 shows an increasing trend in the monthly projected discharge of River Astore under different climatic scenarios. Under the RCP8.5 scenario the discharge increased significantly (p ≤ 0.05) from 196 m³.s⁻¹ to 220 m³.s^−1, whereas a non-significant increase in river discharge was observed starting from 124 m³.s⁻¹ with the maximum value of 151 m³.s⁻¹ under RCP4.5 towards the end of the century.

Mean annual discharge increased by 127.65 m³.s⁻¹ to 136 m³.s⁻¹ during near future scenario (2010–2040); whereas the decreasing trend was observed during the middle future (2041–2070) from 150 m³.s⁻¹ to 145.8 m³.s⁻¹ and the discharge situation remained same at the end of 21^st century under RCP4.5. Whereas under RCP8.5, mean annual discharge increased by 125 m³.s⁻¹ to 150 m³.s⁻¹ during (2010–2040); 160 m³.s⁻¹ to 180 m³.s⁻¹ during middle future (2041–2070) and 210 m³.s⁻¹ to 230 m³.s⁻¹ during far future (2071–2100) projection (Figure 12 and Figure 13).

3.5. Pattern of Hydroclimatic Shift in Astore Basin

Throughout the world, the pattern of the hydroclimatic shift has been reported by many [87,88]. However, its magnitude varies spatially on a regional scale. Temporal changes in streamflow in the Astore basin shows a pattern of seasonal shifts showing a substantial decrease in flow during March–July (Figure 14). On the other hand, it shows a slight increase in discharge during autumn and winter seasons. Historically the Astore basin had a large rainfall-runoff and snowmelt components; however, the reduction in precipitation during spring and the monsoon along with an increase in evapotranspiration will most likely cause a decrease in runoff (Figure 14). The mean projection for RCP8.5 for the period 2071–2100 also shows the same intra-annual pattern for the near future but with more intensity. Many studies [5,84,85] has projected the receding of the glaciated areas in the UIB, which will significantly reduce the share of glacial meltwater that will ultimately cause a reduction in the discharge during the summer season. Also, the decline in the glacier meltwater, our analysis shows a reduction in the monsoon precipitation along with a highly significant (p ≥ 0.01) increase in evapotranspiration in the RCP4.5 and RCP8.5 climate change scenarios. The earlier onset of snow and glacier melt will cause more increase in Astore river flow as indicated by discharge results. This change in hydrological pattern might have strong implication for the cropping pattern and onset of the sowing seasons for different crops in the downstream agricultural areas, especially for Punjab known as the breadbasket of Pakistan.

The results of future forecasted flows are in line with studies conducted by another researcher in sub-basins of UIB [86]. Initially, the summer river runoff will be increased by the mid 21st century and will ultimately decrease at the end of the 21st century. Whereas, for the spring season it will be increased throughout the years. Overall the Astore river flow will be increased, and it is in accordance with the flow of UIB as reported by [87]. These studies used the same RCP 4.5 and RCP 85 scenarios for future forecasting of the hydrology of UIB. However, the comparisons of the results of this study with the results of other studies are subject to the use of different climate models and future scenarios. The uncertainties in hydrological forecasting have also been reported by Prudhomme, et al. [88] at UIB level.

Uncertainties in climate and hydrology projections exist due to the scale of application, data quality and scarcity lead to use of same parameter data for entire complex topography, model forcing etc. The uncertainties exist in the parameterization of modeling are challenging. However, the emphasis should be given to the uncertainties in input data and how they simulate the hydrology of basin. Several scientists [8,76,89,90,91] reported that the rain gauge at high altitude could capture only 20%–30% of precipitation due to topographical variation, blowing and drifting of snow, high-pressure wind. Jonas et al., 2009 reported the range of uncertainty is 13% and 30% in observational data Gigerwald valley. To minimize the uncertainty in observational data Zweifel and Sevruk [92] suggested that wind and evapotranspiration data should be considered for the correction of precipitation. Another factor of uncertainty is that the basin-wide extrapolation of precipitation from the scarce gauges will generate an inappropriate outcome and also the correction of precipitation in high glacierized regions may not be possible due to high altitudinal variation and its dependency on climate forcing [93]. So for an accurate assessment of water resources there is a dire need to have high-quality observational data from dense meteorological stations. The uncertainty of hydrological model does not only depend on the observational data but also on the evaluation criteria of model calibration. In this study model calibration was assessed by three evaluation methods d, R² and NSE coefficients that showed the significantly good agreement with the observed data.

For the downscaling and correction of Global Climate Models GCMs the reference APHRODITE data was used. The gridded APHRODITE data also introduced uncertainty because the data was biased corrected based on a lack of knowledge of glacier mass balance [94]. Finger et al., 2012 stated the glacier and snow dominant catchments and downscaling method of GCMs introduce the significant uncertainty in projected climate change scenario. The uncertainty in climate projection and statistical downscaling method LIBC was evaluated using the QUMP (Quantifying Uncertainty in Model Prediction) ensemble data [95].Temperature lapse rate and precipitation gradient factors for statistical downscaling of RCP4.5 and RCP8.5 reduced the uncertainty in resolving the complex topography of Astore basin. The LIBC technique minimizes the error from a range of 1.7–11.3 °C to a range of 1.4–4.0 °C and negative slopes of linear regression model depicted that error has reduced in the projected downscaled data of temperature as compared to raw GCM data.

4. Conclusions

The present study is the first application of the fully distributed hydrological TOPKAPI model in the Astore catchment, which was focused on the climate projection under RCP4.5 and 8.5 to assess the impact on the hydrological regime of UIB. The TOPKAPI model was successfully calibrated and validated in the Astore river basin using an advance Ordered Physics-based Parameter Adjustment Method OPPA. From the results and statistical measures, it can be concluded that the TOPKAPI model performed efficiently for both the calibration and validation periods. The future forecasting of the hydrology of Astore basin showed a non-significant increase in summer flows and an increased flow during other seasons from 2021 to 2050. Whereas from 2071 to 2100 the summer flows were also found to be decreased with a sharp increase in flows during the other seasons. Overall the increasing trend in Astore river discharge and decreasing trend in snow cover was observed, which can be related to the increase in precipitation (rainfall), temperature and evapotranspiration. The results of this study should, therefore, consider and incorporate into water resources management plans and development of water policy for future water storage in proposed Diamer-Basha Dam and other planned reservoirs on the Indus River system.