Climatic trends variability and concerning flow regime of Upper Indus Basin, Jehlum, and Kabul river basins Pakistan

The Indus Basin is referred to as a “water tower” which ensures water storage and supply to sustain environmental and human needs downstream by a balanced combination of precipitation, snow, glaciers, and surface water. The Upper Indus Basin (UIB) combines the high mountain ranges of the Hindukush, Karakoram, and Himalaya (HKH); this unique region is largely controlled by seasonal meltwater associated with snow and glacier melt during the summer months. The present study seeks to evaluate changes in hydrological and meteorological variable data collected through a network of 35 hydrometric and 15 climatic stations, respectively, across the UIB, Jehlum, and Kabul river basins in Pakistan. The Innovative Trend Significance Test (ITST) in combination with the Modified-Mann-Kendall (MMK) test was used for seeking trends, while Sen’s method was applied for the slope determination of detected trends over four periods of differing lengths (T1: 1961–2013; T2: 1971–2013; T3: 1981–2013; and T4: 1991–2013). Significant decreases were observed in the mean summer and distinct months of (June–August) temperature (Tmean) at most of the stations during T1, while significant increases were dominant over the shorter T4. The mean precipitation (Pmean) was observed as significantly negative at ten stations during July; however, positive trends were observed in August and September. For streamflow, significantly upward trends were observed for mean summer, June and July flows (snowmelt dominant) during T1 and T2, within the glacier-fed basins of Hunza, Shigar, and Shyok; in contrast, streamflow (glacier melt dominant) decreased significantly in August and September over the most recent period T4. For snow-fed basins, significant increases were observed in summer mean flows at Indus at Kachura, Gilgit at Gilgit, and Alam Bridge, Astore at Doyian during (T1–T3). In particular, a stronger and more prominent signal of decreasing flows was evident in T4 within the predominantly snow-fed basins. This signal was most apparent in summer mean flows, with a large number of stations featuring significant downward trends in Jehlum and Kabul river basins. The present study concludes that the vulnerability of this region related to water stress is becoming more intense due to significantly increased temperature, reduced precipitation, and decreasing summer flows during T4.


Introduction
The "water tower" refers to the mountains which ensure water availability and storage to bring about human demands in a sustainable environment (Immerzeel et al. 2020). The Indus Basin is classified as an active Water Tower Unit (WTU) in Global Mountain Biodiversity Assessment (GMBA) due to topographical, elevation, and surface roughness which intersect major river basins (Körner et al. 2017). This basin encompasses the Hindukush-Karakoram-Himalaya (HKH) Mountains and is referred to as High Mountain Asia (HMA), which is a significant source of freshwater for millions of people who inhabit the region (Immerzeel et al. 2010;Lutz et al. 2014). The HKH region contains three large river basins, namely Jehlum, Kabul, and the Upper Indus Basin (UIB) (Hasson et al. 2014) of which the Astore, Gilgit, Hunza, Shyok, Swat, and Shigar are sub-basins. Approximately 80% of the summertime (June-September) flow in Pakistan comes from this region (Ali et al. 2009). Annual runoff of snow and glacier meltwater from the UIB provides around 50% of the recharge for major reservoirs (Mangle and Terbela) and for irrigation, which makes it the biggest continuous irrigation system in the world (Immerzeel et al. 2010;Mukhopadhyaya and Khan 2014;Mukhopadhyay and Khan 2015;Lutz et al. 2014Lutz et al. , 2016. However, flows from the region exhibit large seasonal variability which may result in significant offsets between demand and supply (Reggiani and Rientjes 2014;Jain et al. 2007).
The orientation and topography of the Himalayas result in the interception of substantial amounts of moisture during the Indian summer monsoon season (Jun-Sep) when the seasonal heating of the Eurasian continent results in the advection of moisture from the Indian Ocean. In addition to the summer monsoon offshoots, the precipitation regime combined with mid-latitude westerly disturbances dominates the hydrology of UIB. Moisture transport by the mid-latitude westerly disturbances occurs along the southern flank of the Atlantic and Mediterranean storm tracks originating over the Caspian Sea and the Atlantic Ocean (Pang et al. 2014;Hasson et al. 2015;Bengtsson et al. 2006). Fowler and Archer (2005) and Tahir et al. (2011) previously demonstrated that the Karakoram region was subject to westerly circulation, which results in maximum precipitation being experienced during winter, mainly as snow, whereas the Greater Himalayan region is controlled by the impact of the monsoon regime. Consequently, the Hindu-Kush and Karakoram regions demonstrate different climatic controls on precipitation, compared to the Himalayan Range (Fowler and Archer 2005;Tahir et al. 2011).
Due to the different climatic response of the region, the UIB has been the focus of numerous previous studies. For example, Archer and Fowler (2004) suggested significant increases in annual, summer, and winter precipitation over the period . Subsequently, the authors also detected a significant increase in winter temperature and a decrease in summer temperature within the UIB over the same period (Fowler and Archer 2006). Earlier studies found a significant decrease in temperatures during the summer monsoon months, with positive trends detected during the spring months of March, April, and May (MAM) (Sheikh et al. 2009). Khattak et al. (2011) detected positive trends in winter (DJF) temperatures and a negative trend in summer (JJA) during ; however, they found no evidence of any significant seasonal trend in precipitation. Bocchiola and Diolaiuti (2013) observed increasing trends in winter temperatures at Gilgit and a decreasing trend during summer at Bunji. Similar to Khattak et al. (2011), Bocchiola and Diolaiuti (2013) found no evidence of significant trends in precipitation at Chitral, which lies in the North West Karakoram, or in the UIB.
Several authors have assessed the impacts of climate change on glacier area, melting rates, and streamflow, using a variety of techniques including remote-sensing and hydrological modelling over the UIB (e.g., Akhtar et al. 2008;Ali et al. 2009;Immerzeel et al. 2010, Tahir et al. 2011Butt and Bilal 2011;Latif et al. 2019;Muhammad et al. 2019;Latif et al. 2020b). In contrast to other studies in the region, there are negative trends in precipitation and temperature over the UIB and slightly increasing snowcovered area in Gilgit (Hasson et al. 2015;Latif et al. 2016Latif et al. , 2019Latif et al. , 2020a. The cryospheric studies in the HKH region have found general glacier retreat and a reduction in the snow-covered area based on remote-sensing techniques using satellite data for mass balance computation (Brethier et al. 2007;Kääb et al. 2012). According to these studies, the glaciers in the eastern and central HKH region are losing mass and area and exhibit a negative mass balance, consistent with glaciers elsewhere in the Himalaya's and globally (Bolch et al. 2012;Cogley 2012;Kääb et al. 2012). In contrast, glaciers in the Karakoram region are exhibiting positive mass balances (Hewitt 2005). Hewitt (2005) termed the behaviour of both stable and advancing glaciers in the region as the Karakoram Anomaly, subsequently called the Pamir, Western Kun Lun-Karakoram Anomaly Kääb et al. 2012Kääb et al. , 2015Farrinotti et al. 2020). Varying climatic conditions are attributed to the contrasting pattern of precipitation accumulation between the west and east Karakoram ). In the western Karakoram, westerly circulation results in snowfall during the winter months, whereas the Indian summer monsoon controls the precipitation regime over the eastern region (Fowler and Archer 2006;Bookenhagen and Burbank 2010). As the river flows at the outlet of the UIB are snow and meltwater dominated, they are impacted by the condition of the Karakoram glaciers (Mukhopadhyaya and Khan 2014); there is a direct relationship between river flow patterns and the state of the glaciers in the Karakoram Range. The runoff generated from this region is controlled by snow and glacier melt in the high-altitude zone, while lower-altitude precipitation is responsible for the river runoff (Archer 2003;Ali and De Boer 2007).
The Indus River in the UIB, including its tributaries which are divided into the sub-basins Astore, Gilgit, Hunza, Shigar, and Shyok, is measured at ten stations. Due to the diversity of hydrological regimes linked with snow and glaciers, these sub-basins can be discriminated based on their correlation with climatic variables (Hasson et al. 2015). Previous studies assessing climate change impacts on the UIB (Archer 2003;Fowler and Archer 2006) classified the major basins into the snow-or glacier-dominated regions depending upon the snow, glacier, and runoff contribution to flow resulting from preceding low-altitude winter precipitation. Based on such classification, the Astore and Gilgit basins are classified as snow-dominant basins, while Shigar and Shyok are classified as glacier-dominant basins. According to the Randolph Glacier Inventory (RGI. 4) (Pfeffer et al. 2014), Shigar basin features the highest glacier extent, at 30%, followed by the Shyok basin as 24%; in the UIB, the Astore has 14% glacier coverage. The calculated glacier extents of the Gilgit and Hunza River basins are 9.2% and 28%, respectively (Tahir et al. 2011;Latif et al. 2019Latif et al. , 2020a. Previous studies carried out in the UIB proposed reduced glacial flow contribution while increased annual flow (Sharif et al. 2013) in combination with drying spring and summer cooling (Hasson et al. 2015;Latif et al. 2016Latif et al. , 2020b. However, some hydrological studies concluded increased river flows at the Gilgit River attributed to increased glacier melt contribution to runoff (Latif et al. 2020a). In the present study, the Modified Mann-Kendall (MMK) statistical test in combination with Innovative Trend Statistical Test (ITST) is employed on the available climatic stations (15) and hydrological gauges (35) covering study region. The trend analysis involves determining of significance and direction of a particular trend. Statistical significance is generally assessed using the MK test, while Sen's slope method is employed to determine the direction of a trend. Time series typically display serial autocorrelation which can compromise the trend detection; consequently, several authors suggested removing such effects using pre-whitening and variance correction (e.g., Von Storch and Navara 1999;Yue et al. 2002;Khaliq et al. 2009;Rivard and Vigneault 2009).
Recently, the practice of null hypothesis significance testing (NHST) in hydrometeorological series has come under increasing scrutiny, with the recognition that proposing a hypothesis for a trend based on either external data or physical reasoning and consequently finding significance in a trend from the data themselves is circular reasoning. This was explored in depth by ) who showed that: (a) trend null hypothesis tests are not devised for exploratory analysis of hydrological data; (b) adjusting procedures accounting for correlation in trend tests are insufficient or flawed; and, (c) inductive trend tests cannot provide information on the non-stationary of a process without a priori deductive arguments.
Specifically, Serinaldi and Kilsby (2016) highlighted comprehensive theoretical flaws and ineffectiveness of the use of Trend free Pre-whitening (TFPW) method developed by Yue et al. (2002). They proposed a revised method, by introducing a corrected version called TFPWcu based on their findings estimated by an improved two-stage bias correction of q estimates for autoregressive (AR) (1) signals.
We have carried out a comprehensive trend variability analysis, using the (MMK) method, to overcome TFPW flaws and examine the individual hydrologic response of each basin and magnitude of climatic change over the UIB. This investigation focused on various sub-basins of the UIB by estimating the effect of temperature and precipitation on streamflow, which is highly sensitive to inconsistencies in precipitation and temperature in the form of direct (precipitation) and delayed (snowmelt) runoff. According to the main period of record, during May, the streamflow initiates which gains the maximum level (peak) in July in all sub-basins of the Indus River. This is the fact due to the increased temperature in spring (MAM); the freezing level shifts upward over the winter snow which starts melting and comes to an end in July indicating the end of snowmelt season. The recession flow starts in July in snow-dominant basins (Astore and Gilgit) due to very small glacial extent around 12% and 7%, respectively, according to recent estimates (Yu et al. 2013;Hasson et al. 2014;Latif et al. 2020b). The timing and magnitude of September flows rely on the glacial extent; for the basin having more than 20% coverage, the recession flow starts in early September (Forsyth et al. 2010). Therefore, the present study is mainly restricted to average summer (June-August) and individual summer months for snow-fed basins (Astore, Gilgit, Jehlum). September was included due to its significance in late ablation period for glacier melt in glacier-fed basins (Hunza, Shigar, and Shyok).The objectives of the present study are given as follows: 1. To assess the changes in summer temperature and precipitation trends over the period 1961-2013. 2. To assess the impact of changing climatic trend changes on summer runoff within the glacier and snow-fed basins of the UIB.
These objectives were addressed by examining current trend variability in temperature, precipitation, and streamflow over the UIB, and Jehlum and Kabul river basins. While several previous investigations have been conducted on the UIB (e.g., Archer and Fowler 2008;Khattak et al. 2011;Sharif et al. 2013;Minora et al. 2013;Hasson et al. 2015, Latif et al. 2016Yaseen et al. 2020;Latif et al. 2020a, b), these studies were solely based on conventional Man-Kendall (MK) test providing limited coverage of the available hydrological stations restricted to the UIB excluding Jehlum and Kabul river basins. Moreover, inhomogeneity in the climatic data was also highlighted by Archer and Fowler (2004), specifically in winter minimum temperatures using a double mass curve test. However, this method was only able to indicate the presence of an inhomogeneity rather filling the indicated gap in the available climatic data. Besides, Río et al. (2013) and Forsythe et al. (2014) reported that the records at Gilgit Gupis, Astore, and Skardu, over the period 1952-2009, were homogeneous; however, we identified breaks in the available data at the same stations. The present study expands on previous work, through the inclusion of additional stations, including for the Jehlum, Kabul, and Indus river basin, and undertakes a gap -filling approach for the available data. We employ a state of the art method, i.e., Multiple Analysis of Series (MASH) (Szentimrey 2008), to address the data quality control for making our analysis a reliable trend investigation.
Moreover, to address issues regarding hydrometeorological trend analysis, we used an (ITST) method developed by Şen (2014) which can overcome the issues raised in NHST and TFPW (e.g., Serineldi 2018; Poppick et al. 2016, Clarke 2010. ITST is a non-parametric-based test which is free from any restrictive assumption and follows its application which is rather simple with the concept of sub-series comparisons that are extracted from the main time series  such as T 1 , T 2 , T 3 , and T 4 for the present study.

Study area and datasets
The present study includes the stations within UIB covering western higher and lower regions (e.g., Astore, Hunza, Shigar, Shyok, and Gilgit basins (sub-basins of UIB)), including the Jehlum and Kabul river basins as shown in Fig. 1. The Karakoram Mountains are geographically varied while the Himalayas are considered the third pole due to the large quantities of snow and ice storage, outside the polar region (Farinotti et al. 2020;Li et al. 2018). The UIB is located in the range 33°40′ to 37°12′ N latitude, 70°30′ to 77°30′ E longitude.
We used a digital elevation model (DEM) to delineate the catchment boundary of stations upstream of Khairabad in the Attock district, at the confluence of the Kabul with the Indus River (Fig. 1). The study area lies within the borders of Pakistan (254-8570 m a.s.l.) and excluded the catchments in China and India; however, the stream flows in Pakistan depend on rivers originating in India. The source of the Indus River lies in the Northern Himalayas at the Kaillas Parbat, Tibet (altitude 5500 m), and is joined by few major and minor rivers as it flows through Pakistan.
The streamflow gauging network in the UIB is operated by the Water and Power Development Authority-Surface Water Hydrology Project (WAPDA-SWHP). The UIB contains three foremost basins: the Indus, Jhelum, and Kabul. These basinshave21sub-basins,andstreamflowgaugesareinstalled over the entire UIB at different locations (Fig. 2). The river discharge calculations are based on two steps: calculation of gauge height and measured discharge. The streamflow data from 1961-2013 was obtained from the WAPDA-SWHP. We included the stations having long-term data only after getting excluded stations with less than 20 years of record. This resulted in a total of fifteen climatic stations and 35 hydrological stations being selected for analysis as shown in Table 1. A continuous record of river stage fluctuations is obtained either from periodic readings of a staff gauge height or wire weight gauge or from an automatic water-stage recorder (SWHP Publication 71). Physiognomies regarding each gauging site are given in Table 2. Previously, Mukhopadhyay et al. (2014) and Hasson et al. (2015) used various method for gap-filling techniques to measure the Shigar discontinued discharge, as this gauge  Woo and Thorne (2003) for estimating flow in an ungauged glacierized basin having an area of less than 15,000 km 2 for the remaining period of Shigar discharge as q=Fa -G , where q is specific discharge in m 3 s -1 km -2 , a is basin area in km 2 , F and G are regressed coefficients which are calculated using the flows of gauged rivers.
The data collection procedure for climatic variables by the Pakistan Meteorological Department (PMD) has been comprehensively discussed elsewhere (e.g., Latif et al. 2016Latif et al. , 2020b. The network of hydrological stations for recording and collecting data at river gauging stations is organized by Water and Power Development Authority (WAPDA). Moreover, the Snow and Ice Hydrology Project (SIHP) is a  core project of WAPDA for integrating the hydrological and meteorological data at high-altitude region ultimately making this data available to the end-user (Hasson et al. 2015;Latif et al. 2016). The aforementioned agencies are state-level government organizations who are responsible for the installation and maintenance of the instrumentation, data procurement, and subsequent data broadcasting, following the standard operating procedures (SOPs) compiled by World Meteorological Organization (WMO) (Latif et al. 2016).

Runoff generation in UIB
The mountain ranges of the Hindukush-Karakoram-Himalaya (HKH) region encompassing the Tibetan Plateau (TP) feature a complex high to the low land hydrologic system, including a range of water resources. The significance of the mountains to the total annual flow of the major rivers of Asia, and the sources of runoff within individual basins, varies throughout the region (Yu et al., 2014). The Indus River originates north of the Himalayas on the TP, with headwater tributaries in China (TP), India, Pakistan, and Afghanistan. The main stem of the river flows within Ladakh (Jammu and Kashmir) entering into the northern area of Pakistan Gilgit-Baltistan (GB), flowing between the western Himalayas and Karakoram mountains (Yu et al., 2014). The total estimated Indus river basin inflow which enters from China to India is 181.6 km 3 , while the total estimated inflow entering from Afghanistan to Pakistan in the Indus basin is 21.5 km 3 , out of which 15.5 km 3 is contributed by Kabul River. Interestingly, out of 15.5 km 3 , 10 km 3 is solely contributed by Kunar River which enters in Afghanistan and ultimately flows back to Pakistan. Some small tributaries of Afghanistan (Panjsir, Gomal, Margoa. Shamal, and Kuram) also contribute 6 km 3 to Indus. According to Indus Water Treaty (IWT), western rivers (Jehlum and Chenab) and eastern rivers (Ravi, Beas, and Sutlej) bring mean annual inflow into Pakistan at the rates of 170.2 km 3 and 11.1 km 3 , respectively. However, eastern river inflow is reserved only for India according to the treaty (FAO 2011). The streamflow volume increases along this reach of the river, with contributions from tributaries (gauged) entering the main river from the catchments in the Karakoram Mountains (Shyok, Shigar, Hunza, Gilgit, and Western Himalayas (Astore River) (Hewitt and Young, 1993) as well as ungauged basins from the northern slope of Western Himalayas (Byrne, 2006). Glacier runoff contributes around 24.18 km 3 to the total annual flow of UIB. The Karakoram, Himalayas, and Hindukush contribute around 18% of the total annual flow of 135.6 km 3 over the headwaters of the Indus River at the rates of 14.1%, 2.3%, and 3.2%, respectively; however, the remaining 82% is likely generated from the melting of the winter snowpack (Yu et al., 2014). Immediately north of Mt. Nanga, the westernmost of the high peaks of the Himalayas, the river Indus, turns to the south and flows into the delta of Arabian Sea at Karachi (Sindh Province).

Materials and methods
Here, we focused on a study period of 53 years, from 1961 to 2013, which are then divided in four sub-intervals, namely T 1 (1961-2013), T 2 (1971-2013), T 3 (1971-2013), and T 4 (1981-2013). We analyzed trends of Q mean , P mean , and T mean for summer and determined the number of stations that exhibited positive/negative trends during each data interval. The methods used for the trend detection of the climatic trends and the homogeneity analysis of climate data were discussed comprehensively elsewhere (e.g., Latif et al 2019Latif et al , 2020a. The methodology flowchart for the present study is given in Fig. 3.

Mann-Kendal test
This test has been previously used by various studies (Salmi et al. 2002;Mavromatis and Stathis 2011;Tabari et al. 2011Tabari et al. , 2012.The MK test statistic S is calculated as follows: where x j denotes the sequential data, and n is the length of data (years) for a time series. The mean annual figures during year j and k represent x j and x k , respectively.
Mann-Kendall statistic Z mk , is calculated as follows: Z mk represents positive and negative values representing an increasing and decreasing trends, respectively, at a particular site. After getting confirmed regarding trends, a two-tailed test was used to detect existence at α level of significance either upward or downward pattern persists. In case of the fixed number of Z mk > Z 1−a/2 during required significance stage, the H o was rejected.

Sen's method
We calculated variation during the individual data periods using Sen's method (Sen 1968), based on the following formula: The slope of these N values of Q is figured out using Sen's method. N values of Q i are designated from least to biggest as follows: 3.3 Trends identification and serial correlation effect Kulkarni and Von Storch (1995) previously proposed a method to consider the effect of autocorrelation in time series analysis. This technique was subsequently improved upon by a number of authors (e.g., Zhang et al. 2000;Yue and Wang 2004) and recently applied by Aziz and Burn (2006), Novotny and Stefan (2007), Kumar et al. (2009), andOguntunde et al. (2011). The method calculates the trend and lag-1 auto-correlation of a time series following an iteration until achieving the precise values (Hasson et al. 2015). A number of steps are involved in the pre-whitening of data as discussed in detail by Latif et al. (2016). However, some flaws were mentioned by Serinaldi and Kilsby (2016), particularly regarding hydrological data analysis comprehensive theoretical defects and incompetence of the use of TFPW method developed by Yue et al. (2002). For this purpose, we used MMK test suggested by Hamed (2008) for the removal of serial correlation. The autocorrelation coefficient is computed as where n represents number of terms, x represents lag where "n" is the number of terms; k is the lag; xt is the observed value and is the mean of the first "n-k' terms, and is the mean of the last "n-k" terms. Some authors Cunderlink and Burn (2004) used autocorrelation at lag 1 (ρ1) which is sufficient for hydrological data for testing the hypothesis of serial independence. The null hypothesis (H0) is ρ1 = 0, which means that the data set is serially independent; and alternate hypothesis (H1) is ρ1 > 0, which indicates existence of an autocorrelation. The significance of the hypothesis is tested by Student's t distribution at "n-2" degrees of freedom. The relation is given as In the case of │t │ ≥ t a/2 , the null hypothesis is rejected, and variance is rectified by a factor given as follows where "n" refers to the total number of observations; n*s is an effective number of observations for autocorrelation, and ρk is the autocorrelation function.
Now Eq. 3 will be written as

Innovative Trend Significance Test
The Innovative Trend Significance Test (ITST) method (Şen 2014) is a non-parametric test and consequently does not have any restrictive assumptions. This method is applicable to any time series, with or without serial autocorrelation. The ITST method enables the selection of sub-temporal mid periods for mutual comparison of series which ultimately creates trend on objective and quantitative manners. The approach is implemented as follows: Innovative trend identification methodology follows a linear trend function between independent time variable (t1) and dependent variable (y) by an equation given as where a and b represent intercept on y-axis. For the deterministic methods, parameters are assessed by few alternatives given as follows: 1. The parameters are calculated by eliminating method in which independent (t 1 , t 2 ) and dependent variables (y 1, y 2, ) are known as two points. These two points are substituted in Eq. 8. 2. The slope is calculated by b = (y 2 − y 1 )/(t 2 − t 1 ) which is substituated in Eq. 8. After substitution, only one parameter remains unknown which is then calculated by inserting coordinates of either one of the given point in Eq. 8 as a = y 1− b t1 or a = y 2− b t2 . 3. In innovative trend, any time series is divided into twohalves by sorting each time series in ascending order; finally, two-haves are plotted against each other Şen (2012, 2014). For the construction line, dependent variable sequence values (y 1, y 2, y 3, …. . y n, ) are used which is called "Data line." The brief explanation is given as follows: a. The time series is split into two equal sub-series. b. Each sub-series is arranged in ascending order of rainfall magnitudes. c. The sub-series with earlier observations is plotted as abscissa and that with recent observations is plotted as ordinate, thus making a scatter plot. d. A 1:1 line, i.e., 45°slope, is plotted. e. The point which lies above the 1:1 line exhibits an increasing trend, point lying below the 1:1 line exhibits a decreasing trend. f. The point lying on the 1:1 line does not display any trend.
Using this method, the qualitative trend characteristics of the time series were obtained. The detailed explained in (Şen 2017a, b).

Time series inhomogeneity analysis
The non-climatic issues cause inhomogeneity in time series data. There are multiple reasons which may result in inhomogeneity including recording gauge relocation, and different procedural and observational changes by making discontinuity in time series data which usually compromises accuracy (Aguilar et al. 2003;Costa et al. 2008). These authors stressed to eliminate the inhomogeneity or to determine the introducing inaccuracy which leads to data breaks. These authors anticipated that the non-climatic issues in field, i.e., relocation of detecting gauges and apparatus, setting up improved gauge (station/sensors), altering the observation methods and codes sensor calibration/changing, vary according to local condition.

Standard Normal Homogeneity Test
1. Alexandersson (1986) developed a statistic T (k) for the comparison of the first (k) years of the record with that of the last (n-k) years: 2. When a break is introduced at the year K, then T (k) attains a maximum value near the year k = K. The test statistic T 0 is calculated as: Jaruskova (1994) proposed a similar test for homogeneity analysis using the relationship between test statistic T (n) and T 0 .
3.5.2 Homogenization of climatic record The present study employs the Multiple Analysis of Series for Homogenization (MASH) method proposed by Szentimrey (2008). MASH uses a relative tactic for detection and correction as well as inhomogeneity within time series record. This particular method can be used in a region where metadata is unavailable for a station selected for subsequent trend analysis (Latif et al. 2020b). The MASH enables the direct use of metadata rather taking given values. (Condon et al. 2009).

Field significance
We analyzed the field significance of trends according to the number of flow gauges in each region as discussed in Tables 1 and 2. In order to eliminate the effect of cross/spatial correlation of a station network on the field significance of a particular region, Douglas et al. (2000) proposed a bootstrapping method. This method has previously been used in the analysis of hydrometeorological time series analyses (e.g., Wilks 2006; Guerreiro et al. 2014;Hasson et al. 2015). The approach is given as follows: This method preserves the spatial correlation within a station network but eliminates its influence on testing the field significance of a trend based on the MK statistic S.

Results and discussion
4.1 Temporal and spatial distribution of seasonal streamflow (Q mean ) On behalf of linear regression analysis, Fowler and Archer (2006) allocated the Karakoram region into high-, middle-, and low-altitude stations based on glaciered area, runoff dependence, and precipitation. According to these authors within the high-altitude region (Hunza, Shigar, and Shyok), annual/summer runoff is attributed to mainly summer temperature for glacier melt, followed by middle-(Astore, Gilgit,) and low-altitude gauging stations (Khan Khawer and Siran) in which runoff is solely controlled by winter precipitation (solid) and winter (liquid)/monsoon, respectively. This region can be further divided on such division into the glacierdominant Karakoram Range (Hunza, Shigar, and Shyok) and snow-dominant western Himalaya and Hindukush (Astore, Chitral, Indus, Kunhar, and Swat) (Sharif et al. 2013;Hasson et al. 2015). We have also employed the same division for UIB including Jehlum and Kabul river basin to explicate the results of each stream gauge according to region.

Hunza, Shigar, and Shyok sub-basins of UIB (glacier fed)
We used specific discharge to estimate the discharge variation based on per unit area in glacier-fed basins. The river flows of Hunza River at Daniyor Bridge, Shigar River at Shigar, and Shyok River at Yugo. We observed that the headwater catchments of highly glacierized basins generate high flows in mean summer months, August and September, as shown in Fig. 4. However, high flows during July result mostly from the melting of winter snow cover over the high-altitude regions, while the second peak, in August, exhibits glacial melt associated with the permanent ice sheets and glaciers (Mukhopadhyaya and Khan 2014). Some recent studies (Armstrong et al. 2019) differentiated the sources of meltwater contribution in terms of snow on land (SOL), snow on ice (SOI), and exposed glacier ice (EGI) within Indus Basin. They further explained the melt sources from (SOI+EGI) dominates above 5000-6000-m elevation band in the whole Indus Basin where peak appears in July and August.
We noted a significant downward trend in summer (JJA) runoff in at Daniyor Bridge during the first three-time series (e.g., T 1 , T 2 , T 3 ), at the rates of 4%, 5%, and 8%; however, Shigar discharge increased by 28% during T 1 (Fig. 5). In contrast to the downward trend at Daniyor Bridge for T 1 ,  T 2 , and T 3 , we found significant increases in summer flows of 11% and 7% at Daniyor bridge and Yugo, respectively, during T 4 . These results are only in part consistent with those of increasing long-term trends within the whole of Karakoram (Hasson et al. 2015), because we noted positive trends only in summer flows.
Significant positive trends in discharge were observed during individual summer months (July-August) within all-time series. We observed a significant downward trend in June flow at Daniyor Bridge but an upward trend at Shigar, of 46% and 26%, respectively, during T 1 (Fig. 6). A negative trend at Daniyor bridge is in partial agreement with the findings of Hasson et al. (2015) and Sharif et al. (2013) who found a decreasing discharge at Daniyor Bridge and within the whole of Karakoram region during the  and (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) periods, respectively. We found significant positive trends in July flows at Daniyor bridge during T 3 and T 4 (Fig. 7). These significant positive trends at Daniyor bridge are inconsistent with the negative trends found here during June and July by Sharif et al. (2013). They found an insignificant negative trend during 1966-1995, but significant positive trends during the other months except for April and May. Our results are also inconsistent with the downward July flows within the Karakoram Range found by Hasson et al. (2015). Melting of winter snow represents the main constituent during J u n e a n d J u l y f l o w s w i t h i n S h i g a r w a t e r s h e d (Mukhopadhyaya and Khan 2014). The rising flows are consistent with the positive snow cover trends during winter (accumulation period) in Gilgit basin (Latif et al. 2020a) and the reduced snow cover during the ablation period in Hunza basin (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) as discussed by Hasson et al. (2014). Similarly, we observed consistently significant negative August flows at Hunza River at the rates of 7%, 16%, and 12% during T 1 , T 2 , and T 3 respectively. However, positive trends in August flows at Shigar were observed at the rates of 7% and 21% during T 1 and T 2 , respectively (Fig. 8). Increasing flow at Shigar during June and August is consistent with the previous findings at this station during these 2 months (Mukhupadey et al. 2014). Furthermore, significant decreasing flows at of Hunza River within the Hunza basin are consistent with the summer cooling and increasing snow cover in Hunza and Gilgit basins (Minora et al. 2013;Hasson et al. 2015;Latif et al. 2020a, b).
We observed that September flows significantly increased at Daniyor Bridge, of 19%, during T 2 (Fig. 9). We also noted a positive trend in Shigar flows but that was statistically insignificant. During T 4 , significantly negative flows at the rate of 33% were observed at Daniyor Bridge. The decreasing trends in September flow usually represent glacial melt, which might be attributed to less melting of glacier indicating positive basin storage suggested by recent studies in the Karakoram (Hasson et al. 2015;Sharif et al. 2013). The anomalous behaviour of the Karakoram glaciers is attributed to unique and localized behaviour of weather patterns which makes these glaciers cold and dry even in summer (Kapnick et al. 2014). Similarly, a recent study ) based on the differential Global Positioning System (dGPS) observations in the western Karakoram (Burche Glacier) suggested a slight increase in the debris cover of the glacier attributing to less melting rate correlating with a debris layer thickness (Mayer et al. 2006). However, only decreasing flow can never be enough to establish any linkage between glacier melting rates as Armstrong et al. (2019) suggested only 3% contribution from EGI in Indus Basin which is surprisingly way less as compared to the earlier studies and unable to represent any changes in glacier mass balance. Furthermore, the timing, magnitude, and contribution from the permafrost and debriscovered glaciers particularly in the Karakoram region are still lacking regional-scale calculation which may help to reveal further meltwater sources in a more refined way. The glacier properties in terms of debris cover-thickness, subglacial and surging behaviour of the Karakoram glaciers requires further validation to understand the evolution of anomaly (Farrinotti et al. 2020). Karakoram glaciers exhibit close-to-balance mass budgets which evident the existence of Karakoram Anomaly, however, these signals are more aligned to the regional characteristics rather than an anomalous behaviour (Hewitt 2014). Keeping in view this fact, the most recent field-based glacier study (Muhammad et al. 2020) proposed that a thin debris-layer negligibly impacts the melting rates in the Karakoram glaciers which also reflects another contrasting signal towards the melting as compared to the glaciers where it acts significantly.
The four summer months (JJAS) have underlying importance for identifying the temporal trends because during this season meltwater generates from the glaciers and snowpacks which contain seasonal or perennial snow cover (Mukhopadhyay et al. 2014). We observed a similar pattern of trends in the temporal distribution of stations using ITST, as shown in Fig. 10. The highest number of stations exhibiting negative flows is revealed in summer and September during T 2 and T 4 , respectively. Some authors (Hasson et al. 2015)  associated reduced August flows to less melting of glacial mass in August. Some studies proposed glacier shrinkage in eastern and central Karakoram ultimately exhibiting a negative mass balance (Bolch et al. 2012), endorsing the global trends of glacier shrinkage. On the other hand, positive mass balance and advancing glaciers have also been reported during the last two decade in the Karakoram (Farrinotti et al. 2020). They reviewed conflicting signals that are linked with summer temperature negative trends upward snow cover. Some studies also endorsed the similar stability of glaciers in the Karakoram region, taking into account the glacial mass balance and hydroclimatic trends (e.g., Minora et al. 2013;Sharif et al. 2013;Hasson et al. 2015;Minora et al. 2016). The present study's findings partially endorse the results of the aforementioned research, verifying the negative flows in the Karakoram region: we observed upward trends in June and July flows at Daniyor bridge (western Karakoram), but decreasing flows in August and September indicating lower melting rates of glaciers during ablation. Some potential reasons are associated with the reduced net energy for melting of snow and ice on the Karakoram glaciers such as increased snowfall in the accumulation zone, increased cloud cover and higher snow surface albedo (Farrinotti et al. 2020).

Astore, Gilgit, and Swat sub-basins of UIB (snow fed)
The Gilgit and Astore basins represent Hindukush and western Himalayan ranges in the UIB. The glacial extent of these basins exhibits the least areal coverage like that of Hunza, Shigar, and Shyok basins. The glacier and snow-covered area in Gilgit and Astore basins is 9.2% and 14% establishing 4% and 3% of the entire UIB cryospheric coverage, respectively (Hasson et al. 2014;Pfeffer et al. 2014;Latif et al. 2020b). Despite representing a relatively small percentage of the total area, the glacier coverage contribution to snowmelt runoff should be taken into account to calculate the more reliable runoff estimates (Latif et al. 2020a, b). Therefore, we have classified some flow gauges based on glacial and snow (mixed) flow contribution in the main annual flow, where snow is dominant but glacial flow contributes to some extent, particularly in UIB, according to this classification, Indus at Kharmong and Kachura, Gilgit River at Gilgit and Alam Bridge, Astore River at Doyian falls in UIB, whereas Kunhar River at Naran in Jehlum Basin and Kabul River at Kalam in Kabul River Basin where runoff completely dependent on snowmelt and winter precipitation. The dominant flow in June and July flows is attributed to snowmelt in Astore and Gilgit basins as shown in Fig. 4.
During summer, significant positive flows at Kachura were observed at the rates of 48%, 27%, and 18% during the first three T 1 , T 2 , and T 3 data series, respectively; however, flows decreased by 2% during T 4 (Fig. 5). Similarly, increased flows were observed at Gilgit at the rates of 4% and 6% during T 1 and T 2 . Similar consistently positive trends were also observed in Gilgit River at Alam Bridge at the rates of 10%, 5%, and 18%. Significantly positive flows at Astore River at Doyian during first three data series decreased during T 4 . Similarly, we observed significant positive June flows at Kachura during T 1 , T 2 , and T 3 (Fig. 6). Kharmong flows decreased by 18% during T 1 , but none of the flows was significant in the other three data series. We found a similar pattern of decreased flows in Kunhar river at Naran at the rates of 26%, 42%, 13%, and 8% during (T 1 -T 4 ). Similarly, decreased flows were observed in Kabul River at Kalam at the rates of 27% and 44% during T 1 and T 4 , respectively.
We observed significant decreasing flows during June at Alam Bridge, Kalam, and Khairabad during (T 1 -T 4 ) as shown in Fig. 6. Significant positive July flows at Bunji and Doyian were revealed during T 1 . We found Kachura and Chitral exhibited significant downward flows during T 1 ; similarly, significant negative flows at Kachura were revealed during T 3 (Fig. 7). Similarly, August flows at Kachura increased during T 1 and T 2 as shown in Fig. 8. The seasonal snow at an elevation less than 3500 m starts melting in May causing high flows. The second pulse generates during June due to withdrawing of 0°C isotherm (the altitude at which 0°C temperature maintains for a selected time, i.e., day/month/year) (Mukhopadhyaya and Khan 2014) over the high-altitudinal zones. We observed significantly decreasing flows during August at Alam bridge, Khairabad, and Kharmong during T 2 and T 3 . Similarly, significant decreasing flows were exhibited at Kalam during all respective data series (T 1 -T 4 ).
September flows significantly decreased at Kalam and Khairabad but increased at Kachura during T 1 as shown in Fig. 9. During T 2 and T 4 , Besham Qila, Chakdara, Doyian, Gilgit, Kachura, Kalam, Kharmong, Naran, and Shatial Bridge exhibited significant increasing flows. Some authors observed the August peak instead of recession flow at Besham Qila attributing to glacial melt (Yu et al. 2013). According to our results, this peak has been extended to September due to increased glacial melt at high altitude particularly observed at Astore, Gilgit, and Indus river in Hindukush (Latif et al. 2020b) and Himalayan region. We also noted significant positive flows at Chakdara during T 2 and T 3 , but a significant negative trend in flows was observed during T 4 .
Particularly, T 4 features stronger and prominent signals of a downward trend in flows during 1991-2013. These findings agreed to the findings of Farhan et al. (2014) and Hasson et al. (2015), who also found negative discharge during the last two decades including the Indus River at Kachura. But we found significantly decreased flow over Doyian on the Astore River during the fourth data series, while they reported increased flows over the Astore river during 1995-2012. Only four rivers exhibited increased flows during T 1 : the Indus River at Kachura; Gilgit at Gilgit; Poonch at Kotli; and Brandu at Daggar. The Indus and Kabul rivers show significantly increased flow at Bunji, Besham Qila, Chitral, and Nowshera respectively. This increased discharge might be associated with the observed warming over the UIB during T 3 , though Jehlum basin exhibits a similar pattern of decreased significant discharge.
The highest number of significant and decreasing flows was found during September at 15 out of 35 stations during T 4 . These stations included the Indus, Kabul, and Jehlum river basins, with the highest decrease in discharge over the Indus at Besham Qila and Kharmong, Chakdara at Kabul and Kohala at Jehlum river basin. We suggest the significant drop in snow-fed catchments of the Indus is due to the combined effect of less snowmelt in summer and a significant decrease in precipitation over the UIB (Latif et al. 2019). Similarly, a significant decrease in discharge at Besham Qila (mixed-fed) shows the effect of higher summer cooling. The river flows during the four summer months significantly decreased during T 1 . The streamflows during these months have continued to fall during all data series over the Indus and Jehlum river basins. We noted significant downward trends mostly in the Kabul river basin during June, while positive trends are found only at Astore and Bunji during T 4 . These findings prove the early shifting of recession flows which will result in a further reduction in summer flows. We find spatial variation in the river flows during individual summer months within all data series. Hasson et al. (2015) differentiated snow-fed and glacier-fed regions with the help of a hydrograph based on runoff generation time and suggested that the Indus at Kharmong, Gilgit at Gilgit, and the Astore River at Doyian are basically snow-fed basins showing their peak runoff in July while the rest of the basin is classified as glacier-fed with a peak in August. We found a consistent increase in August discharge for the Indus River at Kachura and Shigar during T 1 and T 2 . These findings are similar to those of Khattak et al. (2011), Mukhopadhyay et al. (2014 and Hasson et al. (2015) who also suggested positive discharge trends at Kachura and Shigar in June and August. Mukhopadhyay et al. (2014) reported an apparent increase in August flow with glacial mass which melts in August, while a rise in September flows occurs due to the glacial melt and monsoonal snowfall over high altitude (> 3500 m). Over the 1961-2013 period, we noted a consistent decrease in discharge from most of the hydrometric stations, which is higher in magnitude and statistically significant.

Jehlum, Siran, Kunhar, and Khan Khwar sub-basin of UIB (foothill catchments)
The low-altitude gauging stations such as Bara, Siran, and Brandu mostly relied on the seasonal contribution of precipitation particularly attributed to the spring season (Archer 2003). Similarly, some authors (Sharif et al. 2013) have divided these foothill catchments based on seasonal snowmelt and monsoon rainfall. We have examined more foothill stations, with the longest and earliest records of flow covering almost the entire Indus, Jehlum, and Kabul basins downstream terrain. Trend evaluation from the past two decades  revealed strong evidence of negative trends of summer flows with high magnitude of decreased flow percentage: 21 gauging stations showed significantly decreased river flows, while only one gauge exhibited as significantly increased discharge as shown in (Fig. 5).
We find a similar pattern of a significant downward trend of flows during individual summer months within all data series. In June, most of the hydrometric stations showed significant downward flows during T 1 as shown in Fig. 6. Seven stations, Chirah, Chahan, Domel, Dhok Pathan, Gurriala, Karora, and Khairabad, exhibited decreasing flows. During T 2 , only Muzaffarabad and Thal exhibited significant downward trends, while Daggar showed a significant upward flow. During T 3 , we noted downward trends in flows at D a g g a r , D h o k P a t h a n , a n d P a l o t e . C h a h a n , Muzaffarabad, Karora, and Thal exhibited significant upward flows. We observed a similar consistent feature of significant negative trends during T 4 at Chahan, Domel, Karora, Kohala, Massan, and Nowshera.
Most of the gauges in the Jehlum basin showed negative discharge in July during T 1 and T 2 ; however, significant upward flows were observed in Kabul River at Chakadara and Jhansi Post. The similar significantly decreased flow was observed during T 4 as shown in Fig. 7.
We noted significant downward trends in August flows at most of the gauges during T 1 and T 2 , except Daggar and Kotli exhibited significantly increased flows as shown in Fig. 8. During T 3 , the number of gauging stations exhibiting negative trends in flows was less as compared to the first two series. Only four gauging stations showed significant downward trends, at Chahan, Domel, Karora, and Palote during T 4 .
Trend analysis for September revealed a pattern of negative trends similar to that expressed in the previous summer months during T 1 as shown in Fig. 9. Chirah, Gurriala, Jhansi Post, Karora, Khairabad, and Palote exhibited significant downward flows. During T 2 Chirah, Chahan, and Gurriala exhibited significant downward trends. During T 3 , Chahan, Chirah, Dhok Pathan, and Gurriala followed a consistent pattern of negative flows, while Jhansi Post showed significantly positive flows. During T 4 , nine stations showed a significant drop in discharge: Azad Pattan, Chakdara, Dhok Pathan, Jhansi Post, Kotli, Kohala, Massan, Palote. Mostly decreasing trends are attributed to weaker monsoonal cycle in terms of drying at low-altitude stations during 1995-2012 (Hasson et al. 2015).

Temporal and spatial distribution seasonal precipitation (P mean )
The inhomogeneous stations and change points in precipitation during 1961-2003 were detected and homogenized using (MASH) (Latif et al. 2020a) as discussed in Table 3. We noted a mixed response of upward and downward trends during the distinct months of summer (Jun-Sep) at different stations installed within the UIB. In June, Dir and Kohat exhibited positive significant precipitation during T 1 as shown in Fig.  6 We found a significantly positive precipitation trend at Dir, Kakul, and Skardu during T 2 . During T 3 Cherat, Chilas, Chitral, and Dir showed significantly positive precipitation while Bunji, Drosh, and Skardu exhibited significant negative precipitation trends. During T 4 , we noted significant upward precipitation at Astore, Bunji, Dir, Kakul, and Skardu while Cherat, Chitral, Drosh, and Gupis showed a significant negative trend in precipitation.
In July, significant downward and upward trends in precipitation were found at an almost equal number of stations during T 1 and T 2 as shown in Fig. 7. Similar dominant significantly negative trends were observed during T 3 and T 4 , except Saidu Sharif and Skardu which exhibited significant upward precipitation. We observed the similar pattern of precipitation trends in the temporal distribution of stations using ITST as shown in Fig. 11(a). In July, T 3 and T 4 revealed the highest number of stations exhibiting downward trends in precipitation.
Interestingly, we observed dominating significantly upward precipitation trends in August during T 1 and T 2 as shown in Fig. 8. Only Dir, Drosh, Kakul, Astore, Chitral, Dir, Drosh, and Saidu Sharif showed significant downward trends. However, the effect of the significantly positive trend in precipitation reduced during T 3 and T 4 , and only Dir and Kakul showed significant upward precipitation in T 4 .
September witnessed a significant upward precipitation trend during T 1 at Astore, Bunji, Chitral, Gupis, and Saidu Sharif while Cherat, Dir, Drosh, and Kakul exhibited significant downward trends as shown in Fig. 9. During T 2 , Astore Bunji, Dir, Drosh, Gupis, and Kakul revealed significant upward precipitation. Dir, Kakul, and Saidu Sharif exhibited significant downward precipitation trends. During T 4 , Astore, Bunji, Gupis, Kakul, and Skardu revealed significant upward precipitation; only Dir and Saidu Sharif showed significant downward precipitation. Our result of upward trends of Astore precipitation during the last two decades endorses results observed by Minora et al. (2013), indicating an increased number of wet days in JAS during 1980-2009. Previous studies reported indefinite and insignificant precipitation over UIB (Khattak et al. 2011;Bocchiola and Diolaiuti 2013) and decreasing annual precipitation (Latif et al. 2016), while upward trends in summer, winter, and autumn but drying in spring at low-altitudinal stations (Hasson et al. 2015). Several authors (Ridley et al. 2013;Cannon et al. 2015;Madhura et al. 2015) have reported consistency between observed winter precipitation and incursions of westerly disturbances. Due to westerlies regime under climate change,  and the number of decreasing wet days, a weakening and northward transfer of the rainstorm trajectory causes increased winter precipitation and decreased spring precipitation (Bengtsson et al. 2006;Hasson et al. 2015). We noted a similar continuous significant drying pattern during the summer (JJAS) when most stations show significant decreasing trends. T 4 revealed an unambiguous drying pattern by exhibiting negative trends (significant) at 8 stations. Astore, Bunji, Chilas, Chitral, Cherat, Dir, and Drosh remain dried during summer although we found a consistent drop in summer precipitation at Chilas, at rates of 8, 10, and 13 mm/decade respectively during T 1 , T 2 , and T 3 . Similar results were found by Hasson et al. (2015) for low-altitude stations in the case of summer precipitation during 1995-2012. However, the decreasing precipitation we found for Chilas was inconsistent with their results since they reported increased monsoonal effect at high-altitude stations including Chilas. We observe a much higher magnitude of wetness in September, the last month of the monsoon period.

Temporal and spatial distribution of summer mean temperature (T mean )
The inhomogeneous stations and change points in mean temperature during 1961-2013 were detected and homogenized using (MASH) (Latif et al. 2020b). Summer temperature significantly decreased at most of the stations during the first three data series but increased during T 4 as shown in Fig. 5. In June, all stations exhibited negative significant trends: however, a significantly positive trend at Saidu Sharif was observed during T 1 as shown in Fig. 6. Similar dominantnegative trends were observed during T 2 , except for Astore, Chitral, Risalpur, and Saidu Sharif which showed significant upward trends. However, we noted dominant significant trends at Bunji, Chilas, Chitral, Gupis, Kakul, and Risalpur during T 3 . During T 4 , six stations exhibited equally significant positive and negative trends. During T 1 in July, 14 out 15 stations revealed downward (significant) trends. Similarly, during T 2 , 10 stations showed significantly decreasing trends except for Chitral, Skardu, Peshawer, and Kohat which exhibited upward (significant) trends as shown in Fig. 7. The cooling pattern was consistent during T 3 by exhibiting in 10 stations the similar negative significant trends except Astore, Chitral, and Kohat. However, during T 4 , the number of stations showing significantly positive trends increased to 7; as compared to the first three series, only six stations revealed significant downward trends. We observed the dominance of warming pattern in summer (JJA) temperature and individual months in temporal distribution of stations using ITST as shown in Fig. 11(b). The number of stations exhibiting significant warming trends in June, July, August, and September were increased during T 4 as compared to other series.
August revealed negative significant trends at most of the stations, while Dir, Kakul, Kohat, and Saidu Sharif exhibited positive significant trends during T 1 as shown in Fig. 8. Similar negative trends were observed during T 2 except Chitral, Saidu Sharif, and Skardu which exhibited significant positive trends. The negative significant trends were consistent during T 3 at eleven stations; only Chitral exhibited a significant positive trend. The similar negative trends were dominant at most of the stations except Chitral, Dir, Gilgit, and Peshawer during T 4 .
September revealed significant downward trends at 12 stations during T 1 except for Chitral and Dir which showed upward (significant) trends as shown in Fig. 9. In T 2 , 8 stations exhibited significant downward trends; only Astore, Chitral, Peshawer, and Skardu revealed significant upward trends. During T 3 , twelve stations exhibited significant downward trends; only Chitral and Kakul showed significant upward trends. However, during T 4 , significant warming trends were observed at ten stations except for Cherat, Dir, Peshawer, and Skardu which revealed significant downward trends.
Spring cooling has increased during the last two decades . We find the same cooling effect on the spatial scale, with Bunji, Cherat, Gupis, and Risalpur experiencing continuous negative trends during all data series. Our findings of cooling at Gupis and Bunji strongly agree with those of Sheikh et al. (2009). They also found a fall in the mean annual temperature at Bunji, Gilgit, and Gupis, but we find significant warming over Gilgit during T 4 . Our trends in T mean suggest significant summer cooling at most stations during the T 1 and T 2 , while during T 4 fewer stations exhibited significant decreasing trends. Comparing the results of trend analysis with Archer (2005, 2006) who investigated trends at low-altitude stations during 1961-1999/2000 and highaltitude stations (Hasson et al. 2015(Hasson et al. ) during 1961(Hasson et al. -2012(Hasson et al. and 1995(Hasson et al. -2012 results are consistent with the summer cooling effect at all stations including Astore, Gilgit, and Skardu, whereas Bunji exhibited a significant upward trend which is inconsistent with their findings. We also found this summer cooling effect decreased during the fourth data series, with an almost equal number of stations presenting increasing significant trends. Hasson et al. (2015) related this warming outcome at low-altitude stations during 1995-2013 to the spring drying having less cloud cover, ultimately upward the number of dry days for the westerly precipitation. Another factor may be the relatively higher magnitude of drying in spring as observed in the present study. According to Hewitt (2005), cloud cover increases with elevation and is much greater for middle and upper ablation zones. We observed significant warming particularly at low-altitude stations (Cherat, Gupis and Kohat) during 1991-2013 that might be associated with the decreased cloud cover. Some authors strongly agree about this summer warming, such as Fowler and Archer (2006), Liu and Chen (2000), and Shrestha and Devkota (2010), who selected the UIB, the Tibetan Plateau, and Nepal, respectively, to conduct their studies. We noted a robust cooling pattern during individual summer months that is restricted to the first three data series. All stations exhibited significant negative trends during T 1 , July and September. Such observations endorse results of Hasson et al. (2015) as we also noted the stronger cooling effect restricted to the first three data series. During the fourth data series, eight stations exhibited significant upward trends. This cooling phenomenon seems to be associated with the stronger dry pattern during July (less cloud cover).
We addressed streamflow, precipitation, and temperature series across the region for all selected stations. For the purpose of making this analysis more reliable and undoubtful, we have also carried out field significance test during last two decades, i.e., 1991-2013 (T 4 ). We noted the field significance based on positive and negative trends for P mean , T mean , and Q mean for all stations. We observed the dominance of negative trends (stream flows) for all variables in Indus, Kabul, and Jehlum river basins. Particularly, negative field significance was observed during summer (JJA) and September, which is completely consistent with decreasing field significance of negative trends in temperature in summer during T 4 . However, trends in precipitation were not consistent with the field significance results as most of the stations revealed positive trends during summer and individual summer months which were not exactly as field significance negative trends in precipitation during T 4 . We observed field-significant decreasing river flows over the glaciated regions of the UIB (Shigar and Shyok) during summer . Hasson et al. (2015) reported similar decreasing discharge during 1995-2012 in July to that in September which is cooler, associating decreasing trend with less summer melting ultimately resulting in positive basin storage as stagnant glaciers. Our results validate previous water resource studies carried out in this region which have suggested significant decreasing trends in hydrometeorological variables in the glacier-fed regions of the UIB Sharif et al. (2013), attributed to the decreased glacier and snowmelt suggesting a positive basin storage . The present study supports the notion of negative trends in UIB runoff and concludes that runoff from UIB rivers is decreasing. Decreased glacier melt during summer appears to be the basic reason, in combination with negative precipitation and temperature trends at high altitude (Forsythe et al. 2017;Latif et al. 2020b).

Conclusions
Using state-of-the-art techniques and comprehensive analysis based on a long-term available record covering four data series, we have revealed a spatial as well as the temporal pattern of prevailing hydroclimatic trends in the UIB. The present study enables the identification of trends using MMK test and ITST, based on a strategic method of field significance for various stations (hydroclimatic gauges) used in the present study. Our results were verified by ITST after getting applied the conventional methods of trend identification. We observed most of the rivers exhibited decreased runoff during the fourth data series over the glacier-fed and snow-fed basins within Karakoram and the Indus, Jehlum, and Kabul river basins. This river flow comprised seasonal snowmelt, glacial melt, and tributary flow collected by a network of stream gauges installed upstream/downstream of Terbella reservoir. Our analysis suggests a positive trend of flows in June and July within the glacial catchment of Hunza at Daniyor bridge attributing to the melting of winter snowfall, while decreases in August and September flows associated with less glacial melt.
We observed a significant decline in P mean and T mean during 1961-2013, including vigorous signs of negative trends seasonally, and observed a consistent drop in the amount of precipitation and river discharge during all data series in summer. Lowerelevation stations like Cherat, Kakul, Risalpur, and Saidu Sharif also exhibited significant negative pattern. Autumn is the driest season during 1961-2013, possibly due to weaker monsoonal precipitation cycle or fluctuating of seasonal precipitation. T mean also indicates significant cooling during the first two series and warming during the last two series. The spatial distribution of trends shows continuous cooling at the northern stations Bunji and Gupis during T 1 . This summer cooling effect decreased gradually, with the least number of stations exhibiting decreasing trends during T 4 , while Cherat, Gupis, Kohat, and Skardu exhibited significant warming during the same period. Our analysis suggests using the ITST method which would be benefited for using hydrological and meteorological time series analysis irrespective of the autocorrelation/serial-correlation effect and any assumption required in trend analysis. Our research suggests that large glaciers across Pakistan are particularly important to meltvolume contributions for the Indus River. Extensive excursions of accessible glaciers are needed to improve the understanding of the relation between climate change and glacier dynamics, with additional and updated observations to include hydrological modeling and glacier mass balance The present study concludes that UIB's vulnerability is growing in terms of significantly increased temperature, reduced precipitation, and eventually decreased river flows during the recent period from 1991 to 2013. The water availability from UIB and other sub-basins will be further compromised in coming decades due to combined effect of climate change and growing population so water awareness should be stressed in downstream areas. The vulnerability can be even stabilized by conservation measurements and efficient use of water according to the requirement.