Assessment of the Climatic Variability of the Kunhar River Basin, Pakistan

: Pakistan is water stressed, and its water resources are vulnerable due to uncertain climatic changes. Uncertainties are inherent when it comes to the modeling of water resources. The predicted ﬂow variation in the Kunhar River Basin was modeled using the statistically decreased high-resolution general circulation model (GCM) as an input for the Hydrologiska Byråns Vatten-balansavdelning (HBV) model to assess the hydrological response of the Kunhar River Basin under prevailing climate changes. The model’s best performance during the calibration and validation stages was obtained with a regular 0.87 and 0.79 Nash–Sutcliffe efﬁciency in the basin, respectively. Under the high-end emission scenario, a 122% increase was expected in evapotranspiration in the rising season of months during the winter period 2059–2079, and such developments were attributed to an immense increase in liquid precipitation and temperature. The model’s output reﬂects a potential for basin stream ﬂow in terms of increasing liquid precipitation up to 182% at the beginning of the monsoon season in the period 2059–2079 in the scenario of high-end emissions. Moreover, the study produced possible uncertainties in high-elevation zones, where the modeling of a catchment can lead to typical snow ablation and accumulation in future projections. This study revealed that the precipitation rate will increase annually, resulting in an increase in the summer stream ﬂow over the basin, though snow is hardly expected to accumulate in the basin’s future climate.


Introduction
Hydrological systems are considered to be excessively critical with respect to their direct effect on economic and environmental development. The hydrological cycle of a basin is greatly influenced by the basin's physical features, climate, and human activity. Many researchers have found that temperature, precipitation, humidity, and wind trends are changing in climatic variables, particularly in precipitation and temperature [1,2]. In recent years, Taylor Eckhardt et al. (2003) have initiated the use of general circulation models (GCMs) to investigate the consequences of future climate change on global water resources [3][4][5][6]. The spatial resolution (200-500 km) for models limits their application to being considered appropriate for the basic level; as such, an excellent spatial resolution is needed for small basins [7]. Researchers have developed complex, quantitative, and computationally efficient techniques for routing basin scales and GCM resolution. With a high resolution of 5-50 km, the regional climate model (RCM) uses the GCM's efficiency while providing high-resolution quality and extensive data on the basin level. Quantitative equations for studying the basin scale and GCM (precipitation and temperature) were developed with  Figure 1). The Kunhar River flows through the southern plateau of the Greater Himalayas. It comes directly from the Kaghan Valley and through Lulusar Lake. It passes through Bata Kundi, Jalkhand, Kaghan, Naran, Balakot, Kawai, and Gari Habibullah, and it meets the Jhelum River at Rara [42]. As far as the water of Kunhar is concerned, it is rich in algal flora and thus harbors a great diversity of aquatic life. The drainage area possessed by the Kunhar River is 2535 km 2 , and the elevation ranges from 600 to 5000 m. The Kunhar is regarded as one of the transboundary area's most significant tributaries [43]. Figure 2 demonstrates the more important topography, such as the contour lines, gradient, delineated sub-basins, and aspects obtained with DEM. This figure shows that the sloping intensity of the Kunhar River ranges from 0" to 76" km. The plateau along the Kunhar River is moderately sloped. It is crucial, however, to note that the majority of the basin is moderately sloped (>8 • and 29 • ) [44]. The large variety of vegetation includes subtropical grassland forests, temperate grassland forests, agricultural land, high mountain glaciers, and snow [45,46]. This challenging topography was categorized into seven key classifications to investigate the Kunhar River Basin's surfaces. Most of the area of the watershed is occupied by snow, agriculture, and forests. This knowledge was obtained from land coverage data from the Kunhar River Basin with a 1-km resolution. agriculture, and forests. This knowledge was obtained from land coverage data from the Kunhar River Basin with a 1-km resolution.

Data
For the current study, climate model statistics (highest temperature, lowest temperature, and rainfall) were also acquired [47]. The SDSM (statistical downscaling model), a system composed of various systematic explanatory variables and meteorological data, was used with the worldwide weather patterns from HadCM3 to find the decline in precipitation and temperature over the period of 2011-2099 [48,49].
While observing the historical data from three climate stations (average precipitation and temperature), this study gathered data from the Pakistan Meteorological Department (PMD) along with data from the Water and Power Development Authority (WAPDA) for two hydrometric stations between 2000 and 2016. Table 1 summarizes the necessary and appropriate geographic descriptions of the hydrometeorological stations.
To ensure data quality, this study identified outlying values (+/−2 standard deviations). In addition, the data were checked to determine whether or not it was possible to proceed over time and space. The Kunhar station was feasible in terms of both time and space. To forecast the same flow in the basin at 0.5°, BCC-CSM2-MR, CAMS-CSM1-0, MPI-ESMI-2-HR, and HadGEM2AO were used with reliable spatial and time resolutions. Intensive geospatial input data were required in order to extract the watershed HBV model. Topography, land use, soil, weather, and hydrological data were collected from different sources/agencies ( Table 2) and used in this model. Digital elevation models (DEMs), land-use/land-related maps, soil maps, environmental conditions, and hydrological details were the main inputs for the HBV model. A 30 × 30 m DEM (digital eleva-

Data
For the current study, climate model statistics (highest temperature, lowest temperature, and rainfall) were also acquired [47]. The SDSM (statistical downscaling model), a system composed of various systematic explanatory variables and meteorological data, was used with the worldwide weather patterns from HadCM3 to find the decline in precipitation and temperature over the period of 2011-2099 [48,49].
While observing the historical data from three climate stations (average precipitation and temperature), this study gathered data from the Pakistan Meteorological Department (PMD) along with data from the Water and Power Development Authority (WAPDA) for two hydrometric stations between 2000 and 2016. Table 1 summarizes the necessary and appropriate geographic descriptions of the hydrometeorological stations. To ensure data quality, this study identified outlying values (+/−2 standard deviations). In addition, the data were checked to determine whether or not it was possible to proceed over time and space. The Kunhar station was feasible in terms of both time and space. To forecast the same flow in the basin at 0.5 • , BCC-CSM2-MR, CAMS-CSM1-0, MPI-ESMI-2-HR, and HadGEM2AO were used with reliable spatial and time resolutions.

HBV Model and Input Data
The current edition of the HBV model for this study was obtained from the University of Zurich. The HBV model adopted in this study can be used on a monthly or nearly daily basis [50]. The simulations performed are of good temporal quality on the daily scale. Input parameters, such as precipitation and temperature, were required to drive the model to obtain climatological values, such as potential evapotranspiration rates and monthly temperatures ( Figure 3). The feedback of the precipitation, either in the form of frozen snow or liquid rain, depends on the temperature during every rising season [51]. With this taken into consideration, every precipitation type was subsequently treated with the soil moisture when appropriate and efficient precipitation occurred, and the resulting value was used as the external leakage for calculation [52]. The remaining portion of the liquid precipitation was instrumental in adding to the soil water reservoir, which could have quantitatively evaporated [53]. Below the surface, samples of substances in the water were found (Table 1). At the channel of the Kunhar River Basin, the core output of the HBV model was water flux and had three sections: interflow, external flow, and baseflow [54]. After incorporating the appropriate hydrometeorological data, the HBV model was developed on a daily time scale (Figure 2). Naran and Balakot, were obtained to be used during the calibration and validation of the model and to further model the watershed's hydrology.  Intensive geospatial input data were required in order to extract the watershed HBV model. Topography, land use, soil, weather, and hydrological data were collected from different sources/agencies (Table 2) and used in this model. Digital elevation models (DEMs), land-use/land-related maps, soil maps, environmental conditions, and hydrological details were the main inputs for the HBV model. A 30 × 30 m DEM (digital elevation model) from the USGS (National Elevation Dataset) dataset [55,56] is shown in Figure 3A; it was used to derive the drainage network of the watershed, sub-basins, and elevation maps for categorization. Land-use data at a spatial resolution of 300 m were obtained from the European Space Agency (ESA) database ( Figure 3B). The land-use data consisted of eight major groups: urban areas (11.26%), grasslands (18.3%), irrigated croplands (13.08%), evergreen forests (6.47%), water bodies (14.39%), rangelands (3.95%), deciduous forests (0.49%), and mixed forests (32.06%). These land-use groups were used to define the landuse classes for the HBV model (Table 3). The land-use data from the ESA were sufficiently detailed to assign the appropriate land-use types in the HBV model. The global IPCC (The Intergovernmental Panel on Climate Change) land soil classes (5 km resolution) from the UN Food and Agriculture Organization (FAO) were extracted to create a land soil file with four soil classes ( Figure 3C) for HBV modeling in the Kunhar River. The land map's physical features included the water efficiency, texture, saturated permeability, mass density, organic carbon, and soil albedo. Long-term climate records from four climate stations (Balakot, Muzaffarabad, Astore, and Naran) during the time frame from 2000 to 2016, which included daily precipitation, daily maximum and minimum temperatures, and daily solar radiation with relative humidity, were obtained from the Meteorological Department (PMD). Since 1961, the WAPDA (Power and Development Authority) has monitored the normal streamflows in the watershed. The regular streamflow records (2000-2016) for the two stream-gauging stations, Naran and Balakot, were obtained to be used during the calibration and validation of the model and to further model the watershed's hydrology.

Hydrological Modeling Method
The HBV model incorporates a routine to calculate daily evapotranspiration (PETm) from monthly values. As inputs, the routine needs the long-term monthly mean potential evapotranspiration (PETm) obtained from the Thorn Thwaite method, which estimates evapotranspiration. Although this method for computing potential and actual evapotranspiration is highly accepted and widely used, it is unfeasible to perform these calculations with large data sets, such as long-term monthly temperature average (Tm) and daily mean air temperature (Td). The daily evapotranspiration was calculated by transforming (adjusting) the PETm through the difference between the Td and Tm and a coefficient C. Bergström mentioned that the adjusted potential evapotranspiration is limited to positive values and is not allowed exceeding twice the monthly average [57].

Snowmelt and Snow Accumulation Module
The snowmelt and snow cover were supposed to be dependent upon temperature in this method. The threshold temperature (TT) was perceived as a critical design parameter for the melting and accumulation of snow at the TT. First, the TT was adjusted to −5 degrees (valid hypothesis) and then modified in the range in which the parameter could be changed. Precipitation in the form of snow occurred when a precipitation situation happened under the TT temperature; therefore, rainfall was expected [58]. The precipitation was considered in the water flow whenever the level of the average temperature rose above the TT. However, when the threshold was higher, both snow and rainfall started to contribute to the drain significantly. The snowmelt strength was measured as the amount of water from Equation (2): where Sm is the snowmelt rate in mm/day, DD is in mm/ • C/day, and T1 is the average daily temperature in • C. As described, the DD shows that the rainwater due to snow increased due to a decrease of 3 • C below the snow limit.

Concept of the Ice Module
A systematic daily analysis technique was used to estimate the melting of the ice, which was ultimately related to the glacier's water absorption, and the daily ratio of snow to glacier ice was computed. The percentage of melted ice in the basin changed spatially and temporally, meaning that the meltwater outflow processes were continuously changed [59], as shown by Equation (3): where Qg is the snow runoff, KGmin is the minimum per day, dKG is the runoff coefficient of the glacier, and AG is a 1-mm measurement parameter. SWE is a measure of the snowpack's water, while S is the snowpack's water content.

Efficient Subcritical Rainfall
The rainfall and accumulation of snow that fell into the water runoff basin were usually divided into two sections: the original surface infiltration that participated and then the donation of external discharge. Depending on the amount of soil humidity during the rainy seasons, a surface module was computed with the HBV model. The coefficient of water content (FC) is defined as the minimum soil moisture storage in the subsurface region. A lower soil moisture content with liquid or solid precipitation indicated an increased the amount of rainfall in the flow. When the amount of water flow approached the value of the FC, the permeability decreased and the rainfall's impact on the flow increased [60]. The precipitation was calculated with the absorption of the actual content of the water flow, as shown in Equation (4): where Deff is the effective precipitation in mm, SM is the net soil vapor in mm, FC is the reduction in the amount of groundwater storage in mm, M is the evaluation of daily rainfall in mm, and β is a variable. β indicates the amount of liquid (M + Sm) water that contributes to runoff. In addition, the coefficient of the runoff increases exponentially when the SM exceeds the FC.

Evapotranspiration Module
As an input for measuring the real evapotranspiration, the long-term mean of the monthly evapotranspiration (EPa; m = 1-12) was used. Afterwards, the reformative potential evapotranspiration was determined with respect to the time interval of the specific transmission period by using the observation that the mean daily temperature varied substantially from the actual monthly average temperature, as shown in Equation (5).
If EPa is adjusted as the potential evapotranspiration in mm, then T is the daily temperature standard in • C, Ta is the monthly mean in degrees, and Z is a variable. If EPa is adjusted, the soil vapor and real outcomes of the evapotranspiration are linked through the PWP (Equation (6)), which shows a relationship between water and the actual evapotranspiration.
Here, Eact is the maximum evapotranspiration in mm. Equation (5) clearly shows that once the moisture of the soil is equal to the PWP, the Eact actually changes at a rate opposite to that of the EPa. The PWP describes the evapotranspiration level of a soil vapor, i.e., when the SM is under the PWP, the Eact is beyond the EPa. The amount of evapotranspiration gradually decreases during soil vapor deposition [61], which is a major flaw of the PWP in Equation (6).

Application of Runoff
In this application, the runoff at a watershed outlet is calculated based on water storage. Water storage was used to simulate the flow of water in a near-surface system. The water storage was used to simulate a quick and steady subsurface system from a material point of view. Transport through the basins occurred at a specific constant speed of capillary action of Qperp. Once the permeable measurement in the upper part of the basin was higher than its L T sensitivity (mm], the overspill quickly followed from the upper Q 0 . The runoff responses of multiple outlets were smoother. The elements of monitoring (K 0 , K 1 , K 2 ) responded to centralized storage operations. The original value of K 0 was considered to be less than that of K 1 to produce quick runoff. The third outlet of (Q 2 ) response was faster than that of the second outlet (Q 1 ), and K 2 had a lower value than that of K 1 , as described in Equation (7): where Q 0 is the near-surface flow in mm/day, Q 1 is the interflow in mm/day, Q 2 is the baseflow in mm/day, Qperp is the percolation in mm/day, K 0 is the field-adjacent influx data factor, K 1 is the interflow-reservoir controller for each day, K 2 is the per-day baseflow storage coefficient, Kperc is a percolation storage coefficient for each day, R u is the principal fluid storage, R u is the secondary storage, R u is the secondary fluid storage rank in mm, and A is a combination of the primary and secondary sources Q A = (Q 0 + Q l + Q 2 ), which was achieved for the total simulated runoff Q A .
Finally, for a single-input reservoir, the time-dependent model is determined as a watershed where the specific runoff Q(t) is proportional to the water saved S(t). Thus, we can obtain a (t) variable based on the coefficient of maximum storage according to Darcy's Law (Equation (10)).
Snow cover and water flow, as well as freezing processes, were simulated in the HBV prerendering. In addition, four-dimensional classification in the South, North, West, and East horizontal aspect categories was also improved in order to offer various effects on the snow and ice streams ( Table 4). The flow of ice was transmitted by the regular historical-grade-day rain and snow routines; furthermore, the opposite behavior of the snow was attributed to the reduced surface temperature and the acceleration of the degreeday measurement to remove water [62]. The overall vegetation cover for the three cover areas showed an area of 5682 km 2 at an elevation of about 5320 m.a.s.l.; the highest entire area of 424 km 2 was at approximately 3600 m.a.s.l., and there was prevalent glaciation over approximately 1700 km 2 between 6300 and 5200 m.a.s.l. ( Table 4). The genetic algorithm and Powell (GAP) method was applied to the optimal calibrated values of the HBV parameters in forest land, pasture land, and glacier land, as presented in (Table 4).

Calibration and Validation of the HBV Model in the Kunhar River Basin
The genetic algorithm and Powell (GAP) optimization approach was used for model calibration. The GAP is a stochastic design method and works evolutionarily by selecting and recombining high-performing parameter sets with each other. The GAP algorithm consists of two steps [58,59]. First, optimized parameter sets are generated by an evolutionary mechanism of selection and recombination of a set of initial, randomly selected parameter sets (again within user-defined parameter boundaries). During the second step, parameter sets are fine-tuned using Powell's quadratically convergent method, as described in [60]. The HBV model, as for the calibration method, was used to calibrate the HBV model by the GAP method. With the GAP evolutionary method, the procedure was repeated 100 times to run the model and 1000 times additionally with Powell local optimization [61]. When running the GAP procedure 100 times, 100 different parameter data sets are computed, and then the appropriate set regarding the selected goodness-of-fit criteria is computed. Moreover, the GAP procedure was repeated to run 100 times, but with only 16 parameters (GAP16). In this way, it was tested whether the parameter values were closer to the initial ones rather than by using all 16 parameters (GAP16). A similar fixing of some of the parameters to avoid over-parameterization was done as described in [62]. The model was first calibrated for each of these time periods and then tested using validation periods selected to reflect specific climate conditions (dry, wet, warm, and cold). Each calibration period was validated by three different types of PE input data: precipitation, air temperature, and potential evapotranspiration. The analysis encompassed the entire study period, and the validation periods did not overlap. The procedure is schematically presented in Figure 1. The hydrological model with input data and with known model parameters shown in Table 1 was calculated. The HBV model's performance with the set of parameters was calibrated (2000-2009) and validated (2010-2016) for Gari-Habibullah and Naran, which accurately represent the Kunhar River Basin's hydrological characteristics ( Table 5). The principal origin of the river that flows into the Kunhar watershed is ice, glacier runoff, and rainfall. For this reason, the criteria of the parameters were selected based on these river flow origins. The HBV model frequently ran throughout the calibration and validation for Gari-Habibullah and Naran and worked consistently well between the calibration and validation stages, with Nash-Sutcliffe (NSE) efficiency values of 0.87 and 0.79 around the basin. The relationship between the simulated and observed river flows in the calibration and validation periods became extremely good, with NSE values of 0.87 and 0.79, respectively (Figures 4 and 5). The model reflected excess discharge in 2006 in Gari-Habibullah and high discharge in 2007 during the calibration process, whereas in 2011 and 2015, the discharge was unpredictable and significantly overstated throughout the validation period.        The year 2013 was reported as one of Pakistan's warmest times. The highest average temperature in northern Pakistan was 47.5 °C on 19 July 2013. The temperature might have increased stream flows through excessive ice melting and melting of glaciers, but this was not effectively identified by the model due to undervaluation of the river flows The year 2013 was reported as one of Pakistan's warmest times. The highest average temperature in northern Pakistan was 47.5 • C on 19 July 2013. The temperature might have increased stream flows through excessive ice melting and melting of glaciers, but this was not effectively identified by the model due to undervaluation of the river flows in 2013. By comparison, the years 2011 and 2015 were also mentioned in Pakistan's unpredictable weather history, which is the justification for the accurate assessments; these may be irregularities in the existing meteorological data. Likewise, for the validation period, in Pakistan, 2013 was one of the warmest years, and the temperatures in different regions of the country ranged from 47 to 51 • C. Briefly, the causes behind our overestimations of the projected stream flows may be due to mistakes made with the climate index method in the simulation system, inadequate seasonal rainfall, or extreme precipitation in certain seasons. Similarly, in a past analysis, the seasonal river flows throughout the Kunhar River Basin were uncertain [8]. The hypothesized inadequate seasonal flow in the Kunhar River Basin was focused on environmental impacts relative to several other inadequate catchments in the Upper Indus Basin (UIB), Pakistan, which are dependent on winter rainfall and early spring flow in the Kunhar River Basin. Furthermore, various conditions influence the melting of ice and glaciers, including the landscape, weather, dimensions, elevation, hillsides, and soil types of the region, whether intentionally or unintentionally. After all, the weather scale approach does not accurately describe certain parameters [63]. PBIAS revealed unfavorable results of 0.47 and 14.61, indicating that the calibration and validation cycles were inadequate for the actual model ( Table 5) To evaluate the efficiency coefficient, a statistical approach was performed to determine the hydrological system's predictive control with the Nash-Sutcliffe (NSE) efficiency coefficient [51]. It is defined as where Qobs is observed, Qsim is simulated, and Re f f is the NSE. A reference value of 1 shows an absolutely ideal match. Calibration and validation, which resulted in NSE values of 0.87 and 0.79 around the basin, were highly effective in the HBV model ( Figure 3). After achieving these parameters for the HBV model, the expected river flow changes for the periods of 2017-2037 (near), 2038-2058 (mid), and 2059-2079 (far) were incorporated into the HBV. This was based on four GCMs and four RCP scenarios.

Comparison in Each Time Scale (Daily/Monthly/Seasonal)
At the daily scale, the statistical results show no significant difference between the GCMs rain data, with CC of 0.39/0.48 and AAE values of 14.3/16.4, respectively (Table 6). A difference is found in that the GCMs' values underestimated the actual precipitation, with a PBIAS value of −15.24%, while RCPs overestimated the actual precipitation with a PBIAS of 98.4%. Therefore, the AE value of the GCMs is also much higher than that of the RCPs, with values of 7.3 and 5.8 mm/day, respectively. As predicted, the trend described above is also seen on the monthly scale. The CC values of the RCP data ranged from 0.39 to 0.48, showing a good correlation with the GCMs' data. Moreover, the AE and AAE values of the GCMs are many times higher than those of the RCPs. The errors at the daily scale were eliminated by the aggregation to the monthly scale, causing the CC to become more balanced; however, this does not explain the big difference observed in the evaluation trends between GCMs. The GCMs' precipitation data is always overestimated across the basin, and the largest bias statistic indicator values were recognized with this dataset in evaluations by Mou Tan [65]. Generally, the GCMs' precipitation data are slightly more accurate and agree relatively better than the RCPs' data with observations measured on the monthly scale. As shown in Table 6, the analysis results of the seasonal statistical indicators obtained from the RCPs' data show the largest mean errors, with AE and AAE values that are too large. At the same time, the PBIAS value of the GCMs in the dry season is −39.8%, many times different from the rainy season value of −7.3%. This is related to the very low GCMs' rainfall that occurs in the dry season; the lower rainfall value in the denominator of Equation (7) will cause the PBIAS value to be higher. Due to underestimating rainfall in the dry season, the rainfall in the GCMs' data makes the difference between the two seasons much larger than the observed data. The rainfall ratios between the dry and rainy seasons for GCMs were 17% and 93%, respectively.

Projected Changes in Precipitation
Expected increases in the monthly precipitation in GCMs of the near, mid, and far future for the four RCPs were estimated through a model simulation ( Figure 6). In the warm weather and in early fall, precipitation significantly increases when using CAMS-CSM1-0, showing an increase of up to 58% for return intervals in a deep-term period and 38% percent for the progressing monsoon months in the mid-future period in the RCP 2.9 emission scenario. In the RCP 4.1 emission scenario, BCC-CSM2-MR has an increase of 60% in early monsoon precipitation in the approximated future, although CAMS-CSM1-0 indicates a reduction of approximately 48% in the autumn period in the long term. HadGEM2AO has a forecast for a 61% rise in the beginning of the monsoon cycle for the near future, while CAMS-CSM1-0 predicts a 47% decrease in the advancing monsoon season duration for the upcoming future timeframe in the RCP 6.3 emission scenario. BCC-CSM2-MR estimates an incredible 139% shift in the start date of the monsoon. For the relatively near-future period, a significant 93% decrease is assumed according to the MPI-ESMI-2-HR project in the RCP 8.8 emission scenario. The general opinion in the identification of relevant plans is that there will be an indicative liquid precipitation increase (from the rainy season) and an increase in the intensity of rainfall (from the cold season) in the future, according to the results from the watershed [66].

Projected Changes in Temperature
The predicted periodic temperature changes for four RCPs of the GCMs are shown in Figure 7. The temperatures of the Kunhar River Basin are predicted to rise, on average, until around the end of the decade. According to MPI-ESMI-2-HR, an increase of up to 3.5 °C was recorded in the RCP 2.9 emission scenario during summer-to-winter transition months in the future. According to the MPI-ESMI-2-HR data, a much greater increase of 4.5 °C is expected during the winter-to-summer transition months. In addition, cold-to-warm transition seasons are predicted with MPI-ESMI-2-HR in the RCP 6.3 emission scenario to rise 5.6 °C higher than RCPs 2.9 and 4.1 throughout the winter. Moreover, during the cold-to-warm transition period, an alarming 8.5 °C increase (the maximum increase among all RCPs), as indicated by MPI-ESMI-2-HR, would occur in the future under the RCP 8.8 scenario. All of the simulations of the temperature assume that the period for which the snow stays is likely to reduce the glacial period's metamorphic phase due to the massive variations in all of the emission scenarios across the watershed. In other words, due to the start of mid-summer that is predicted for the basin, a significant glacial reduction is expected [67].

Projected Changes in Temperature
The predicted periodic temperature changes for four RCPs of the GCMs are shown in Figure 7. The temperatures of the Kunhar River Basin are predicted to rise, on average, until around the end of the decade. According to MPI-ESMI-2-HR, an increase of up to 3.5 • C was recorded in the RCP 2.9 emission scenario during summer-to-winter transition months in the future. According to the MPI-ESMI-2-HR data, a much greater increase of 4.5 • C is expected during the winter-to-summer transition months. In addition, cold-to-warm transition seasons are predicted with MPI-ESMI-2-HR in the RCP 6.3 emission scenario to rise 5.6 • C higher than RCPs 2.9 and 4.1 throughout the winter. Moreover, during the cold-to-warm transition period, an alarming 8.5 • C increase (the maximum increase among all RCPs), as indicated by MPI-ESMI-2-HR, would occur in the future under the RCP 8.8 scenario. All of the simulations of the temperature assume that the period for which the snow stays is likely to reduce the glacial period's metamorphic phase due to the massive variations in all of the emission scenarios across the watershed. In other words, due to the start of mid-summer that is predicted for the basin, a significant glacial reduction is expected [67].

Projected Changes in Evapotranspiration
The increases in quarterly evapotranspiration expected were described for four RCPs of the GCMs (Figure 8). For the mid-and long-term winter peaks, the highest enhancement of 46% was reported under the MPI-ESMI-2-HR emission scenario with RCP 2.9. In the highest recorded winter months, a greater increase of up to 64% can be seen for the far-future timeframe, as defined by the MPI-ESMI-2-HR in the RCP 4.1 emission scenario. Furthermore, a significant increase of up to 69% was observed in the coldest months for future periods predicted in the RCP 6.3 emission scenario for MPI-ESMI-2-HR. Furthermore, an impressive 98% change was recorded for the coldest months in the extended period, as shown by MPI-ESMI-2-HR in the RCP 8.8 emission scenario.

Projected Changes in Evapotranspiration
The increases in quarterly evapotranspiration expected were described for four RCPs of the GCMs (Figure 8). For the mid-and long-term winter peaks, the highest enhancement of 46% was reported under the MPI-ESMI-2-HR emission scenario with RCP 2.9. In the highest recorded winter months, a greater increase of up to 64% can be seen for the farfuture timeframe, as defined by the MPI-ESMI-2-HR in the RCP 4.1 emission scenario. Furthermore, a significant increase of up to 69% was observed in the coldest months for future periods predicted in the RCP 6.3 emission scenario for MPI-ESMI-2-HR. Furthermore, an impressive 98% change was recorded for the coldest months in the extended period, as shown by MPI-ESMI-2-HR in the RCP 8.8 emission scenario.

Projected Changes in Stream Flow
The maximum flow variability in the RCP 8.8 emission scenario was calculated ( Figure 9). With BCC-CSM2-MR, monsoon flow is expected to rise by up to 42% in the mid future, while for MPI-ESMI-2-HR, it is projected to decrease by up to 19% in the winter months in the long-term RCP 2.9 emission scenario. The monsoon flow in the next few months showed a rise of up to 63%, and in the extreme winter months, it showed an increase of up to 34% in the middle of the next half according to the IPSL-CM5A-LR state in the RCP 4.1 emission scenario. CAMS-CSM1-0 forecasted a rise of up to 76% in progressive rainy season stream flow over several months, while the RCP 6.3 emission scenario projected a decrease of up to 21% in the peak cold months in the long term. In addition, CAMS-CSM1-0 projected an increase of up to 149% during monsoon initiation periods for the far-future period, and with MPI-ESMI-2-HR, the fall flows are expected to drop to 36% in the RCP 8.8 emission scenario. It was observed that all GCM results are compatible in the direction of flow changes in each emission scenario. Simultaneously, a series of small variations were reported for the river streams during the forecast periods. The effects of the water flow measurements indicate a significant increase in hot season water flows and a complete decline in cold season water flows. This could be due to the

Projected Changes in Stream Flow
The maximum flow variability in the RCP 8.8 emission scenario was calculated ( Figure 9). With BCC-CSM2-MR, monsoon flow is expected to rise by up to 42% in the mid future, while for MPI-ESMI-2-HR, it is projected to decrease by up to 19% in the winter months in the long-term RCP 2.9 emission scenario. The monsoon flow in the next few months showed a rise of up to 63%, and in the extreme winter months, it showed an increase of up to 34% in the middle of the next half according to the IPSL-CM5A-LR state in the RCP 4.1 emission scenario. CAMS-CSM1-0 forecasted a rise of up to 76% in progressive rainy season stream flow over several months, while the RCP 6.3 emission scenario projected a decrease of up to 21% in the peak cold months in the long term. In addition, CAMS-CSM1-0 projected an increase of up to 149% during monsoon initiation periods for the far-future period, and with MPI-ESMI-2-HR, the fall flows are expected to drop to 36% in the RCP 8.8 emission scenario. It was observed that all GCM results are compatible in the direction of flow changes in each emission scenario. Simultaneously, a series of small variations were reported for the river streams during the forecast periods. The effects of the water flow measurements indicate a significant increase in hot season water flows and a complete decline in cold season water flows. This could be due to the essential rise in spring temperatures, which contributed to the monsoon months and significant winter evapotranspiration changes in the GCMs. After all, the measurements of the rapidly increasing flows substantially exceed the amounts of decreasing flows in all estimates that were compatible with the results [68].
Water 2021, 13, x FOR PEER REVIEW 18 of 24 essential rise in spring temperatures, which contributed to the monsoon months and significant winter evapotranspiration changes in the GCMs. After all, the measurements of the rapidly increasing flows substantially exceed the amounts of decreasing flows in all estimates that were compatible with the results [68].

Discussion
In the mountain region of the Kunhar River basin Pakistan, the large amount of precipitation leads to a substantial accumulation of snow during the winter and spring season. The stream flow increases rapidly with the beginning of the snow melting period in March and April. The Kunhar River Basin receives its maximum precipitation in winter and spring, which results in a clear temporal separation in accumulation of the snow and the time of the peak flow. The Kunhar River is characterized by a seasonal cycle of river flow with maximum flow in the summer season. Summer runoff of the upstream

Discussion
In the mountain region of the Kunhar River basin Pakistan, the large amount of precipitation leads to a substantial accumulation of snow during the winter and spring season. The stream flow increases rapidly with the beginning of the snow melting period in March and April. The Kunhar River Basin receives its maximum precipitation in winter and spring, which results in a clear temporal separation in accumulation of the snow and the time of the peak flow. The Kunhar River is characterized by a seasonal cycle of river flow with maximum flow in the summer season. Summer runoff of the upstream catchments of the Southern Asian Rivers is controlled by the melting of snow and glaciers, which, in glacierized catchments, contributes up to 50% of the seasonal flow [69]. The downscaled precipitation and maximum/minimum temperature estimates from the five GCMs were in good agreement with the measurement-based data under the baseline condition. The product of multiple GCMs of the Intersectoral Impact Model Intercomparison Project (ISI-MIP) was used to assess the influences of the climate variation on the hydrological regimes, considered a suitable representative [48,70,71].
The meteorological data from four climate stations and streamflow data from one hydrological station were used in this study. We analyzed the possible consequence of projected climate variation on interseasonal and annual water flow, as well as on snowmelt in the Kunhar River basin. We studied the responses of hydrological processes, such as discharge and extreme and median flows, to climate variation. The temporal impacts of the projected climate change were assessed by coupling a well-calibrated semidistributed hydrological model (HBV) with the results of the five GCMs under two greenhouse gas emission scenarios (RCP 4.1 and RCP 8.8). The effect of topography was corrected by applying the elevation band (temperature lapse rate and precipitation lapse rate) approach in this study, which improved the simulation results. We presented the results of the streamflow calibration and validation for the Kunhar River. Additionally, by using the HBV tool, we studied the sensitivity and prediction uncertainty of the model parameters, which was necessary to evaluate the strength of the calibrated model. The results of the GCMs presented continuous warming over the Kunhar River basin at the annual and seasonal scales in the middle and at the end of this century.
This finding is consistent with the results of previous studies in the South and Central Asian regions, including the Pamir-Alay [72,73], the Tian Shan-Pamir-North Karakoram [74], the Himalayas [75], and the Tibetan Plateau [76]. It is expected that the wet and dry seasons in the basin will become more severe than those in the baseline period. We showed that these phenomena are predicted to be stronger under RCP 8.8 than under RCP 4.1. A possible reason for these alterations is that in the future, more meltwater will be produced in early summer and more snow will be replaced by rain. The typical changes will accelerate the convergence of water flows and raise the flooding frequency and intensity. Folini et al. [77] reported that aerosol emissions in the 20th century might increase, in association with the enormous population and industrialization growth. Similarly, Bollasina et al. [78] confirmed that in Asia, the concentration of atmospheric aerosols has increased steadily.
Xin et al. [79] reported that over China and Central Asia, a rising trend in the concentration of aerosols in the atmosphere could cause a substantial rise in temperature. The mean annual precipitation over the Kunhar River basin is expected to rise in the 2038-2058 future time period under RCP 4.1, as evidenced by the results of two GCMs (CAMS-CSM1-0 and HadGEM2AO), as well as the two GCMs (MPI-ESMI-2-HR and HaDGEM2AO) that indicated a rising trend of precipitation in the second time period (2059-2079) under RCP 8.8. The remaining GCM models that were analyzed in this study showed a decreasing trend of mean annual precipitation in the two future periods under both RCPs. The winter mean precipitation had a lower decreasing trend than that of the summer and fall seasons. The summer mean precipitation exhibited a greater decreasing tendency than the other seasons in the two future periods. Similarly, in the Yellow and Xin River Basins in China [80], in the Middle East [71], and in the westerly dominated region of South Asia, a decreasing trend of summer precipitation was indicated [81]. The increasing or decreasing propensity of winter precipitation varies from model to model. Two GCMs showing a rising trend of winter precipitation and resembling this analysis were reported by Luo et al. [82] for the Heihe River Basin and by Omani et al. [83] for the Pamir-Alay Mountains in Central Asia. However, the patterns of seasonal variations in precipitation for three GCMs presented in the Kunhar River Basin are contrary to those in reported for the Hunza River Basin of the Karakoram Mountains [84] and the Jhelum River Basin of the Himalayan Mountains [74]. Pendergrass et al. [85] reported that the global winter precipitation increased over the second half of the 20th century, and they attributed this to the role of increasing moisture counteracted by weakening circulation. Li et al. [86] pointed out that in Central Asia, at the end of 20th century, there was a persistent decreasing trend of annual precipitation, and Meng et al. [87] confirmed that precipitation might increase in the middle of the 21st century in the south of the Tian Shan Mountains. The multimodel ensemble result projected a decrease in average annual precipitation during the 2038-2058 and 2059-2079 time periods under both RCPs in the Kunhar River basin. The interannual and seasonal scale analyses of the mean precipitation changes presented large uncertainties among the GCMs, as evidenced by both increasing and decreasing tendencies under RCP 4.1 and RCP 8.8 in various future time periods. In the Asian region, these contradictions in our findings may be associated with the rising concentration of anthropogenic absorbing aerosols and the westerlies system [73,74]. The analysis of the simulated Kunhar River flows mostly described an increasing trend of the mean annual streamflow during the 2059-2079 time period under RCP 4.1 and RCP 8.8, as well as a decreasing trend during the 2038-2058 time period under RCP 4.1. Mostly, the increase and a lesser decrease of annual future flow may be attributed to the similar projection of the different GCM models for the total annual precipitation. Similarly, a decreasing/increasing tendency for the different future periods of annual river flow was indicated in an arid alpine catchment in Karakoram [81]. In addition, it is expected that the average monthly peak discharge may shift to earlier in the summer season, from July to June, for the two future time periods under both RCP 4.1 and RCP 8.8 for almost all GCMs, which is mainly due to the slight rise in precipitation in the spring and winter seasons, as well as because of an earlier snowmelt caused by global warming. Similarly, due to earlier snowmelt, Sorg et al. [88] projected the impacts of climate changes on flow seasonality and concluded that less water will be available in the summer months in the Syr Darya River Basin in Central Asia. Olsson et al. [89] confirmed, from trend analysis, this shift in seasonality of flow and predicted a possible decreasing trend of annual flow in the Zarafshan River Basin in Central Asia. A similar shift in the discharge peak (July to June) was pointed out by Liu et al. [90] for the Yarkant River Basin in Central Asia; in contrast, Babur et al. [91] reported that the discharge peak could be delayed (July to August) in the Jhelum River Basin. These discrepancies in detections might be related to various projected climate models for winter and summer seasonal precipitation in Central Asia. We found that the average monthly peak discharge in the Kunhar River Basin indicated a significant decreasing tendency in August and September for the 2038-2058 and 2059-2079 time periods under RCP 4.1 and RCP 8.8 for all five of the GCMs. In this study, most of the GCM model outputs along with the multimodel ensemble results showed that the summer water flow in the Kunhar River is expected to increase at the end of the 21st century under the two studied greenhouse gas emissions scenarios. Increases in summer water flow in the Kunhar River can be ascribed to the rapidly melting snow and ice caused by the continuously increasing air temperature. Meanwhile, it could lead to uncertainty in the predictions due to fluctuations in the reduction of sea ice and the gradual increase in the snow lining. Furthermore, these conditions are predictable if the temperature in all GCMs increases significantly.

Conclusions
This study is the first to refer to the role of temperature in extreme events based on regulations in the Kunhar River Basin. The temperature verification in the Kunhar River Basin shows that RCP scenarios can be representative as ground temperature measurement stations in meteorological and hydrological studies. Temperature events, such as very cold, damaging cold, strong sun, and scorching hot events, affect the rainfall distribution and the inputs to the flow simulations. Moreover, the proposed study is suitable for humid climate conditions in the tropics, such as the climatic conditions in the study area, and can be reliably used in other basins with similar conditions. After incorporating all of the appropriate hydrometeorological data, an HBV model was developed on a daily time scale. First, the model was calibrated and validated using parameters from 2000 to 2016. The predicted flow variation in the Kunhar River Basin was calculated using statistically decreased high-resolution GCM performance as an input for hydrological HBV models under different RCP emission scenarios. The level of maximum rainfall was almost 1650 mm in the wind-retrieval and fall months, with a maximum temperature of 23 • C during the starting months of the rainy season.
The maximum evapotranspiration was almost 4 mm during the starting months of the rainy season, and the highest discharge of flow was about 1032 m 3 /s along the basin; these values show the yearly climatology of the hydrometeorological parameters. The HBV model for stream flow projection included a snow melting and snow accumulation module, ice module, efficient subcritical rainfall, evapotranspiration module, and an application of runoff. This highly effective HBV model was used for calibration and validation, with NSE values of 0.87 and 0.79, respectively. Postprocessing is expected to significantly improve the predictions of liquid precipitation and consistency (winter) in future precipitation slices. The temperature predictions of all of the simulations suggest that the snow residence period is likely to decrease for the metamorphic cycle of glaciation due to significant improvements in all of the basin's emission scenarios. The results also show that the summer streamflow over the basin is expected to significantly increase, as shown by the analysis of the simulations' performance. The advantages and disadvantages of study suggest that local knowledge/information is also very useful in hydrometeorological research to avoid excessive misunderstandings of gridded climate scenarios and can instead be understood as an opportunity to explore the potential of reanalyzing data in terms of their performances that are, as of yet, unproven due to limited, short duration, and heterogeneous observational data. Our tentative studies will be further expanded with other gridded climate scenarios already recognized in Pakistan, and the spatial variations in water balance components and the effects of climate change on flow changes in the Kunhar River Basin will be calculated. This study's findings are essential because emission-based hydrometeorological simulations of the Kunhar River Basin and the multi-GCM output differences across the river's watershed have been insufficiently studied.