Characterization of interannual and seasonal variability of hydro-climatic trends in the Upper Indus Basin

High-resolution seasonal and annual precipitation climatologies for the Upper Indus Basin were developed on the basis of 1995–2017 precipitation normals obtained from four-gridded datasets (APHRODITE, CHIRPS, PERSIANN-CDR, and ERA5) and the quality-controlled high- and mid-elevation station observations. Monthly precipitation is estimated through the anomaly method at the catchment scale, and then, it is compared with the observed discharges over the 1975–2017 period for verification and detection of changes in the hydrological cycle. Running trends and spectral analysis on the precipitation gridded dataset were performed. The Mann–Kendall test was employed to detect the significance of trends whereas the Pettitt test was used to identify change points in precipitation and discharge time series. The results indicate that the bias corrected CHIRPS precipitation, followed by the ERA5, performed better in terms of RMSE, MAE, MAPE, and BIAS against the rain gauge observations. The running trend analysis exhibits a slight increase in annual precipitation, but it shows significant increase in winter precipitation. A runoff coefficient greater than one, especially in the glacierized sub-catchments of Shigar, Shyok, Astore, and Gilgit, indicates that precipitation is likely to be underestimated and glacial melt provides excess runoff volumes in a warming climate. Streamflow variability is found to be pronounced at the seasonal rather than at the annual scale. The annual discharges at Shyok, Gilgit, and Indus at Kachura gauges are slightly significantly increasing. Seasonal discharge analysis reveals more complex regimes, varying in different catchments, and its comparison with precipitation variability favors a deeper understanding of precipitation, snow-, and ice-melt runoff dynamics, addressing the hydroclimatic behavior of the Karakoram region and some weaknesses in the monitoring network at high altitude.


Introduction
Changes in the hydrological cycle at both global and regional scales are a major concerns as these changes occurring either due to climatic or anthropogenic factors can have adverse impacts on society in general. This emphasizes on monitoring these changes on all relevant scales. Precipitation is a key component of the hydrological cycle and perhaps the most difficult one to be monitored, particularly in the glacierized and snow-fed catchments such as the Upper Indus Basin (UIB), where complex Himalayan terrain and harsh environment makes the high-altitude precipitation monitoring further difficult. As a result, the observations are sparse, short, and discontinuous with large gaps, all making it difficult to understand the variability and changes and its subsequent impact on the basin hydrology and water yield (Archer and Fowler 2004). Hence, a combined study of both precipitation and discharge provides a better insight into the variability of the water cycle in the high mountain areas and the cryosphere in general (Poloczanska et al. 2018) and in the UIB in particular.
The UIB spans over the confluence of the Hindukush, Karakoram, and Himalayan (HKH) mountain ranges that constitute one of the largest cryosphere reserves outside the poles (Soncini et al. 2015). The UIB fulfills water demands of rapidly increasing population and feeds one of the world's largest irrigation systems within the Indus Basin. Pakistan's economy to a large extent depends upon the agriculture, and in turn, on the water resources of the UIB, which feeds around 16.5 × 10 6 ha of irrigated land. The UIB up to the Besham Qila inflow gauging station drains an area of 163,528 km 2 and receives 26% of its annual flow from the snowmelt and 44% from the glacial melt (Mukhopadhyay and Khan 2015). The primary usage of the UIB waters is for irrigation and hydropower production as the climate of lower Indus basin is mainly arid and hyper-arid, unlike other South Asian basins, which are characterized by wet regimes of summer monsoon. Immerzeel et al. (2015) marked the UIB as a climatic hotspot region, due to the wide variation in the climate anomalies, as well as due to the significantly rising water demand downstream.
The UIB is characterized by diverse hydro-climatic regimes with contrasting patterns across sub-basins and over time. Archer and Fowler (2004) found rising rates in summer, winter, and annual precipitation. Khattak et al. (2011) found summer cooling and winter warming, and no confirmed changes in precipitation over the 1967-2005. Immerzeel and Bierkens (2012 noted the highest vulnerability of water scarcity conditions for the Indus River Basin among ten basins in Asia. They recognized that significant population increase, groundwater depletion, climate change, snowmelt, and ice melt are key factors that affect the hydrological regimes of the Indus basin. Bocchiola and Diolaiuti (2013) found a slight increase in annual precipitation over northwest Karakoram and Chitral-Hindukush despite, its decrease on the greater Himalayan side. They also observed summer cooling and winter warming, which are more prominent than earlier particularly at Bunji and Gilgit stations. Ali et al. (2015) examined current and future climatic and hydrological changes over the UIB. They showed that northern parts of the UIB experienced a larger increase in temperature and precipitation than southern parts. Projections of future changes show a consistent increase in temperature and precipitation. The rate of increase of river flow is greater in winter compared to summer season. They considered higher river flow possibly due to a larger increase in the air temperature and the consequent enhancement of the melting of the snow and ice cover. Latif et al. (2018) explored both seasonal and annual precipitation trends in the UIB using low and mid altitude stations. The results exhibited significant falling rate of annual precipitation in six stations, while three stations showed a rising rate of precipitation. Overall, the UIB experienced a downward trend in precipitation both spatially and temporally.
The studies mentioned above mainly used valley-based sparse and fragmented low-and mid-altitude stations being operated by the Pakistan Meteorological Department (PMD). The estimates from these stations neither represent the climatology of high altitude areas, nor provide any quantifiable mechanism that draws logical inferences between low-and high-altitude precipitation (Khan and Koch 2018). As a great amount of the UIB streamflow originates from the active hydrological zone in the 2500-5500 m altitude range, data from low altitude stations (even if making up long time series of observations) are not representing reliable hydro-meteorological conditions over the frozen UIB water resources (Hussain et al. 2017).
Moreover, the Indus is a transboundary basin (see Fig. 1) and observed hydro-meteorological data are mostly scattered, discontinuous, and not easily accessible. Hence, it is difficult to assess the spatial and temporal variability in high altitude mountainous regions using a sparse ground-based observation network, as it cannot depict horizontal and vertical precipitation variability effectively (Lutz et al. 2014). Various gridded datasets have been developed based on satellite-based data (Huffman and Bolvin 2013), interpolated observation (Yatagai et al. 2012), and reanalysis data (Baudouin et al. 2020;Li et al. 2020) to handle this issue. Most of past studies concerning the UIB depend upon regional or global-gridded datasets for mass balance and hydro-climatic studies (Baudouin et al. 2020;Dahri et al. 2016;Immerzeel et al. 2009;Iqbal et al. 2019;Krakauer et al. 2019;Lutz et al. 2014;Masood et al. 2020;Minallah and Ivanov 2019;Rizwan et al. 2019;Ullah et al. 2018).
Gridded datasets provide a better solution in terms of temporal and spatial coverage, even if they are affected by likely greater errors, in particular in high altitude areas where large bias and uncertainty may occur, especially in conditions of significant snowfall (Andermann et al. 2011) and in glacierized catchments (Wortmann et al. 2018). Uncertainties and biases in the gridded datasets are usually due to shortcomings of data sources and generation algorithms of these products (Sun et al. 2018). It is also noted that gridded datasets such as CHIRPS and APHRODITE mainly used the World Meteorological Organization's Global Telecommunication System (GTS) gauge data in their production mechanism (Yatagai et al. 2012). The WMO GTS collaborates with PMD for sharing observed meteorological data. In 1995, Water and Power Development Authority (WAPDA) in Pakistan collaborated with the International Development Research Centre, Canada, to install automated weather stations within the Pakistan side of the UIB known as the data collection platforms (DCPs). However, the data of these stations are neither publicly available, nor they are shared with the WMO. Keeping in mind the above issues, the primary objective of this study is to discuss biases and uncertainties in each available gridded dataset (CHIRPS, ERA5, PERSIANN-CDR, and APHRODITE), before and after their corrections performed using mid-altitude PMD and high-altitude WAPDA-DCP stations. The second objective is to compare catchment scale precipitation to discharges to identify possible sources of errors and climatic anomalies, as well as their changes over time and space.
Previous assessments on hydro-climatic trends at the annual or seasonal scale considered specific periods Bolch et al. 2012;Fowler and Archer 2006;Hasson et al. 2017;Janes and Bush 2012;Krakauer et al. 2019;Latif et al. 2018;Masood et al. 2020;Rahman et al. 2018;Sharif et al. 2013;Zaman et al. 2020a, b), but they neither incorporated variations within specific temporal sub-periods, nor they described non-linear dynamics of hydro-climatic variability. In an agriculture-based country, with substantial climatic variation, it is difficult to understand the hydro-climatic phenomenology using conventional linear trend schemes. Building on the need of generating future scenarios of water sustainability, this study shows the results of a running trend analysis that covers the entire dataset period to assess hydroclimatic variability from decadal to interdecadal trends in the UIB at the sub-basin scale. More specifically, this study provides a comprehensive investigation of the detectable links of short-term (sub-decadal and decadal) and long-term (multidecadal) precipitation and runoff variability at the seasonal and the annual scale. In this way, it supports a more detailed overview of precipitation and runoff regimes at the basin and the sub-basin level and improves the past analysis for examining hydro-climatic behavior.
Within this context, a precipitation climatology  for the UIB at the seasonal and the annual scale was built using the anomaly method applied in (Crespi et al. 2018. This method is briefly recalled in the second section of the paper after the description of the study area for the collected precipitation and runoff data collection and of the implemented statistical analyses. In the third section, the hydro-climatic trends in each sub-basin of the UIB are examined, together with biases and uncertainties of the gridded precipitation data compared with the ground observations. A discussion of the results, also in comparison with those coming from the analysis of runoff data, follows in the fourth section.

Description of the study area
The UIB, including part of the western Himalayas, Karakoram, and northern Hindu Kush mountains, lies in the Fig. 1 The location of the Upper Indus River Basin (UIB-light blue in the inset) and of the hydro-meteorological stations used for the analysis geographic domain within 31-37° latitude and 72-82° longitude. The basin hosts 14,085 small and large glaciers, covering an area of 19,338 km 2 (RGI Version 6.0). The drainage area of the Indus Basin is around 163,528 km 2 at the Besham Qila gauging site. The basin is shared between Pakistan, India, and China with around 46% of its area lying within Pakistan's administrative boundaries (Hasson et al. 2017). The Indus River originates from Mount Kailash in western Tibet at an elevation of 5486 m and has an overall length of 3180 km measured at the outlet into the Arabian Sea (Jain et al. 2007). The main stem flows initially through the Ladakh district in Jammu and Kashmir and afterwards it enters northern Pakistan (Gilgit-Baltistan), between the Himalayas and the Karakoram range. The catchment area and, consequently, the discharge of the Indus River become larger in Gilgit-Baltistan (GB) when tributaries, such as Shyok, Shigar, Hunza, and Gilgit Rivers in the Karakoram Mountains and Astore River in the western Himalayas merge with the main river stem. Afterwards, it turns towards south from Nanga Parbat (8126 m asl) and flows through three provinces of Pakistan, i.e., Khyber Pakhtunkhwa (KPK), Punjab, and Sindh, before its confluence to the Arabian sea. Additionally, Chitral, Swat, and Kabul rivers that are originating from the Hindu Kush Mountains also join the mainstream of the Indus river in the KPK province, whereas western Himalayan rivers of Jehlum, Chenab, Ravi, and Sutlej join main the Indus river stem at Punjnad in the Punjab province.
The basin's highest elevation is set by the K2 (Karakoram-2) peak, also known as "Godwin-Austen," the second highest mountain in the world (8611 m asl), whereas the lowest altitude at Besham Qila is 542 m a.s.l. The mean altitude is around 3750 m asl, whereas 35% of the area lies above 5000 m asl. The hypsometric curve of UIB is shown in Fig.S1. The UIB has seven subcatchments, i.e., Gilgit, Hunza, Astore, Shigar, Shyok, Shingo-Zanskar, and Indus Downstream (Mukhopadhyay and Khan 2015). The Shyok and Shigar basins feed the eastern and the central part of the Karakoram. Almost one third of the Shigar basin is covered with glaciers, including the world's largest glaciers and ice masses after polar regions. The hydro-climatic characteristics of each sub-basin are quite different. Summer monsoon and westerlies are dominant sources of annual precipitation in UIB; however, the effect and contribution of both sources vary spatially, as well as temporally (Hasson 2016).
The average annual precipitation measured at different stations within the basin ranges from 156 mm at Gilgit station to around 1514 mm in Pir Chanasi valley. The mean annual discharge at Besham Qila is 2405 m 3 s −1 (Hasson et al. 2017). The Indus Basin receives 70% of its annual flow from June to September with a maximum value in July. October to March are distinguished as low flow months. UIB climate falls into the "cold desert" category (BWK) according to the Köppen-Geiger climate classification, i.e., an area with a little precipitation and a large daily temperature range. The relationship between station elevation and measured precipitation is represented in Fig. 2, which shows that no significant altitudinal trend can be observed, although doubts arise about the reliability of precipitation data when snowfall occurs.

Hydroclimatology of the Upper Indus Basin
In the Upper Indus Basin hydrology, the precipitation regime features annual-round midlatitude western disturbances. Such disturbances sometimes carry a solid form of moisture, mainly during winter and spring (Hewitt 2011). The rate of such unusual solid form of moisture is higher during the positive phase of the North Atlantic oscillation (NAO), when western disturbances affect Afghanistan and Iran due to low heat over the area, resulting in extra moisture input from the Arabian sea (Hasson 2016;Hasson et al. 2014;Syed et al. 2006); Initially, water supply is generated by the melting of snow starting from the mid of March to late June. The extent of water availability mainly depends upon the concurrent temperature and accumulated snow amount (Hasson et al. 2014(Hasson et al. , 2013. Afterwards, snowmelt runoff is merged with glacier melt runoff from late June to late August as a consequence of high air temperatures. The climate of the UIB is classified by winter extra-tropical cyclonic/anticyclonic circulations (westerlies) and South Asia summer monsoon atmosphere circulations. Both winter and summer have a significant impact on the climatic patterns of the UIB (Hewitt 2011). The westerlies enter the UIB through the northwest by the end of November or early in December. Initially, these westerlies are presented in distorted and diffuse states. Afterwards, they interact with the already existing orographic trough with low pressure that allows them to recover their potency and frontal structure. The topographic blocking separates these westerlies into southern and northern sections around the western Tibetan Plateau and Karakorum (Pang et al. 2014). The relationship between topography, local climate, and circulation system determines the net precipitation and distribution pattern in the UIB. The differential heating between land and sea is the main reason for summer precipitation (Dahri et al. 2016). The summer monsoon carries moisture from the Arabian sea that moves along the Indus valley towards the western Himalayas. It also brings moisture from Bay of Bengal and moves northward to the eastern Himalayas, and from Indian ocean to the western Himalayas following the path along the Indus river valley Hasson 2016;Pang et al. 2014).
It is generally believed that the precipitation rate is increasing with elevation up to a certain elevation where a maximum is reached. For the UIB some studies (Immerzeel and Bierkens 2012;Immerzeel et al. 2015;Winiger et al. 2005) found that precipitation is increasing up to a specific elevation, i.e., 5000 m, and shows a downward trend above this elevation. Most of the annual precipitation falls in the winter and spring and originates from the westerlies. Although summer carries occasional rain to trans-Himalayan areas, it accounts for only one third of total annual precipitation (Madhura et al. 2015). Some glaciological studies mentioned a significant increase of precipitation rates of 1500-2000 mm at 5500 m (Soncini et al. 2015). At the high elevation zones, ice is the primary source of hydrological regimes, followed by snow melt, while the contribution of summer monsoon precipitation is small.

Meteorological data
Data availability is a big issue in the HKH especially in the UIB where stations are neither densely nor uniformly distributed. There are primarily two organizations: WAPDA and PMD who are responsible for the collection and management of hydro-metrological data across northern areas of Pakistan. Approximately, hundreds of gauging stations are installed in various places across UIB (Zaman et al. 2020a, b). In this study, meteorological data from 26 stations were selected from PMD, WAPDA, and China Meteorological Data Sharing Network (CMDSN) based on completeness, temporal duration, and homogeneity. Out of these 26 stations, seven are operated and maintained by PMD. These are valley-based stations, which are located within the altitude range of 1200-2200 m asl. The data from these stations were collected for the period 1981-2017. The second meteorological network is maintained by the Snow and Ice Hydrology Project (SIHP) of the WAPDA, which is operating 12 automated weather stations known as data collection platforms (DCPs), located in the elevation range of 1479-4440 m asl and providing observations since 1995. As the Karakoram range hosts the largest snow ice reserves of the UIB, DCP stations operated by the WAPDA are particularly relevant for examining the hydro-meteorological conditions prevailing over the UIB cryosphere (Hasson et al. 2017). In addition, there are also three high altitudes stations operated and maintained by EvK2-CNR (Italian-based organization). Two stations of Askole (3015 m asl.) and Urdukas (3926 m asl) provide observations since 2005, and one located at Concordia (4690 m asl.) provide observations since 2011. All three stations are used to calculate standard meteorological parameters. However, data time series from these three stations feature large data gaps and are affected by uncertainties especially during the winter season due to sensor inefficiency in extreme weather conditions. The detailed information about climatic stations, their elevation, time period, and mean precipitation are given in Table 1.

Streamflow data
Upper Indus Basin is fed by three sources of streamflow, i.e., glacier melt, especially in Shyok, Hunza, and Shigar subbasins, followed by snow melt mainly in the Gilgit and Astore subbasins and rainfall-runoff. The daily streamflow of nine hydrometric stations within the UIB was taken from the Surface Water Hydrology Project (SWHP) of WAPDA, Pakistan from 1973 to 2017, except for Indus at Kharmong and Bunji stations where, data are available from 1983 to 2017 and 1973 to 2013, respectively. Discharge data of Astore at Doyian and Gilgit river at Gilgit are used in this study. Similarly, discharge of Hunza basin is collected at Dainyor station. Table 2 provides specfic information about these streamflow stations and their outflow points.

Gridded observations
In the last decades, a great progress was made in developing analyzed fields of precipitation over regional and global scales providing different gridded climatic products. These products are available at a regional and global scale and are used in hydro-climatic assessment studies. Precipitation products can be divided into four major categories: (1) climatic model reanalysis, (2) satellite estimates, (3) merged satellite and station observations, and (4) rain gauge-based observations (Sun et al. 2018). In this study, we used at least one product from each of these categories in order to check their accuracy for hydro-climatological studies based on precipitation estimates through these datasets. The APH-RODITE is based on station observations, CHIRPS is a combination of satellite and station observations, ERA5 is a reanalysis dataset and PERSIANN-CDR is based on remote sensing using artificial neural network. The detailed information about these products is given in Tables 3 and 4. For instance, APHRODITE (V1101 and V1101EX_R1) is specifically developed for summer in the Asian region with spatial resolution 0.25° × 0.25°; the products are provided by the Data Integration and Analysis System (DIAS, Japan) that is based on an interpolation of 3500 to 8000 gauge observation (Dile and Srinivasan 2014).
The second dataset is CHIRPS which was developed by the Climate Hazards Group (CHG) and the United States Geological Survey (USGS). It is available from 1981 to present with a spatial resolution of 0.05° and temporal resolution at daily and monthly scale. It was developed by combination of ground-based gauge information and cold cloud duration measurement by the synergistic use of satellite infrared radiometers. Passive microwave and GridSat-B1 satellite data were employed to update the PERSIANN algorithm to estimate daily precipitation. It is based on remotely sensed information combined with an artificial neural network (Ashouri et al. 2015). The European Center for Medium-Range Weather Forecasts (ECMWF) launched the new reanalysis product ERA5 data. The analysis is developed using an advanced 4Dvar assimilation scheme at temporal and spatial scales. It is available at 0.25 × 0.25° and computes various atmospheric variables at 139 model levels for 1979-present time period at different temporal scales (Baudouin et al. 2020).

Anomaly method: the interpolation scheme from rain-gauge network to regular grid
The precipitation station observations in the UIB are sparse and do not provide complete temporal and spatial coverage. Therefore, it is inappropriate to develop basin-wide annual and seasonal precipitation climatology based on available observations directly. For this purpose, monthly precipitation records for four gridded datasets (CHIRPS, PERSIANN-CDR, and ERA5 from 1995 to 2017 and for APHRODITE from 1995-2015) were selected, according to data availability. The gridded data were reconstructed over the study area by means of the anomaly method as described by Crespi et al. (2021). To develop the precipitation climatology on the seasonal and annual scale, monthly gridded and observed data were aggregated at the seasonal scale, i.e., winter (DJF), spring (MAM), summer (JJA), autumn (SON), and annual scale. In this scheme, the precipitation   (Saha et al. 2010;Tarek et al. 2020) signal is reconstructed by superimposing the spatial fields of seasonal climatology to the spatio-temporal fields of relative anomaly, i.e., deviations from the reference for a certain season. The ratio or multiplicative correction factor between referenced and gridded precipitation at each ith specific station is computed as below: Here, p test m and p reference,i are gridded and observed precipitation series at the specified time scale over the period of common data availability. It has to be pointed out that the rain gauge observations are here assumed to be "true" reference precipitation values, and the other precipitation products are adapted to them. P test m,I is the seasonal ratio anomaly or correction factor which is then interpolated as described in Crespi et al. (2018). Similarly, the interpolation method is also applied on gridded datasets on the same scale. These two fields are calculated individually, and the season estimates are finally obtained by their product. The same procedure is applied on an annual scale. The anomaly-based climatology helps to develop fine-scale information provided by gridded datasets and incorporates the available records on a large area when the station distribution over the study domain does not provide complete coverage (Brunetti et al. 2012).

Evaluation criteria
To evaluate the data quality of precipitation anomaly against observed precipitation for an overlapped period of 1995-2017 and 1995-2015 at the seasonal and annual scale, four statistics, i.e., bias (BIAS), mean absolute error (MAE), root mean square error (RMSE), and mean absolute percentage error (MAPE) were assessed in this study.  199519951998-2007(Lutz et al. 2014) 1998(Dahri et al. 2016) 1998-2009 199519951998-2007(Lutz et al. 2014) 1971-2004CWC and NRSC, 20142000 17 199519951998(Dahri et al. 2016) 1998-2007(Lutz et al. 2014 where PS i is the value of gridded precipitation estimate for the ith event, PO i is the value of rain gauge observation for the ith event, and N is the number of precipitation events. Bias represents the average difference between the gridded precipitation, after the anomaly method correction, and gauge precipitation, as reference. A negative value of bias implies underestimation, while a positive value indicates an overestimation of observed precipitation. The gridded precipitation datasets were also compared with observed streamflow, and the runoff coefficients were computed as shown in Table 6 in order to assess the ability of runoff data to close the water balance in each sub-basin of UIB and also to verify the reliability of rain gauge-based precipitation assessed at the catchment scale as pointed out, for instance, by ).

Precipitation and runoff trend analysis
The anomaly-derived precipitation and streamflow discharges were used for the trend analysis in sub-basins of UIB at seasonal and annual timescales. The discharge period from 1973 to 2017 was selected for four rivers (Astore, Indus at Besham Qila, Indus at Shyok, Indus at Kachura), 1973-2013 and 1983-2017 were selected for Indus at Bunji and Kharmong and Shigar, respectively, due to the available data. There are two basic tests which are mainly used to detect trend significance, i.e., parametric and nonparametric approaches (Zaman et al. 2015(Zaman et al. , 2016. The trend for this hydro-climatic time series was estimated using robust nonparametric regression techniques, i.e., Mann-Kendall (MK) test (Kendall 1948;Mann 1945) in conjunction with Theil-Sen (Sen 1968;Theil 1950) slope method to determine the trend's slope. Mann-Kendall is a rank-based method that determines the presence of any trend irrespective of the type of sample data distribution and whether such trend is linear or nonlinear (Tabari and Talaee 2011). There are two reasons for using MK test for trend analysis. Firstly, it is not necessary for time series data to be normally distributed. Secondly, it is insensitive to missing values and data outliers and less sensitive to breaks caused by inhomogeneous time series (Bocchiola and Diolaiuti 2013). In a running trend with a moving time window approach, a type of exploratory data analysis similar to that adopted by (Brunetti et al. 2012) is used to calculate and visualise trends for precipitation and discharge over different time windows and access their significance using consecutive years of the datasets as starting point (x-axis in Figs. 4 and 5) and ending points (y-axis in Figs. 4 and 5). As trends in climate change studies are expected to be analyzed after 20 years of monitoring, at least, a minimum duration of 15 years was considered for precipitation trend analysis, due to limited availability of data time series, and 20 years for discharge. However, generally, running trend analysis does not require a fix threshold on the length of time series and the threshold value can be altered according to study objectives, climatic parameters, data availability, and local issues. Such analysis is not only helpful to detect non-linear hydro-climatic trends in UIB over the different temporal scales, but it also facilitates a comparison of these results with other studies which did not show overall climatic fluctuations in the study period.

Change point analysis
The nonparametric Pettitt test is also employed in this study to observe change point in hydro-climatic time series (Palaniswami and Muthiah 2018;Zhang et al. 2015). It supposes a time series X t with t = 1,2….N has a change point at time step T. The values of X t for t = 1,2,…,T have cumulative distribution function (CDF) F 1 (x) and t = T + 1,T + 2,…,N have CDF F 2 (x) and F 1 (x) ≠ F 2 (x). Like other statistical measures, the null hypothesis (H o ) depicts the absence of change point against the alternative hypothesis (H 1 : change point present). Given the random variable k(T) defined as the Pettitt statistics K is written as, And time at which change occurs in time series is determined by p values of the two-tailed Pettitt test is computed by as The null hypothesis (H 0 ) would be rejected for p < α, where α is the significance level. In our study, we keep the significance level at 5%

Climatology, anomalies, and precipitation records
The annual and seasonal precipitation climatology for the study basin is shown in Fig. 3 In summer slightly wetter conditions occur, followed by the monsoon season which receives the major portion of precipitation. Considering the mean seasonal precipitation over the whole study basin, the precipitation rate is 96 mm in winter (DJF), 150 mm in spring (MAM), 220 mm in summer (JJA), and 70 mm in autumn (SON). It is generally believed that the precipitation rate is getting higher with elevation until an orographic optimum. However, the observation of this phenomenon is quite uncertain in UIB where higher precipitation values are measured at lower altitudes or no significant increase with altitude is observed (see Fig. 2). There can be two main possible reasons: the first one is related to the fact that the majority of meteorological stations are located in low lying areas and so there are significant chances of under catch snow precipitation, a common problem in windy and snow-dominated areas like UIB (Petäjä et al. 2016). The second is related to the representation of precipitation with gridded datasets which are often lacking sufficient as well as reliable gauge observations. The values of error statistics (2) to (5) computed with the leave-one-out method over the grid points with gauge observations are presented in Table 5. The bias-corrected CHIRPS precipitation had the best performance at the annual and seasonal scales followed by ERA5 whereas PERSIAN-CDR had the worst performance at both scales.
The summary performance of selected gridded datasets is as follows: CHIRPS datasets is slightly underestimated before correction while ERA5 is highly overestimated before correction. However, after correction, the statistical results of ERA5 are closer to CHIRPS-gridded datasets. The higher values of MAE, BIAS, and RMSE for ERA5 before correction were due to incorporating liquid and frozen water, consisting of rain and snow that fall on the earth's surface. Secondly, reanalysis datasets mostly rely on coupled numerical models, ocean, and atmospheric data and do not depend upon ground-based observations (Copernicus Climate Change Service 2017). The performance of the Aphrodite dataset is found to be inferior to ERA5 especially at the annual scale. A great variation of PERSIAN-CDR with observed values might be associated with bias adjustment, which is based on the GPCC dataset with a coarse resolution of 2.5° resolution (Fallah et al. 2019). The reliability of CHIRPS and ERA5 were also cross-checked by comparing with some previous studies as shown in Table. 4. The results show reliable agreement for developing precipitation climatology using observed and gridded data series extended in space and time in UIB. The discrepancies with other studies are due to varied coverage of study area, time period, and number of meteorological stations applied in the analysis.

Variability and trends of the precipitation in Upper Indus Basin
The 1995-2017 annual and seasonal precipitation records for the study domain were evaluated for short-and long-term trends by using the Theil-Sen slope test (Theil 1950) while the trend significance was evaluated by MK test (Mann 1945). By assuming a confidence interval (C.I) 0.05 (95%), the variability of precipitation on a finer time scale was calculated using running-trend analysis or moving average window approach on the 1995-2017 records. Moreover, change point analysis of precipitation at annual and seasonal time series was also investigated using the Pettitt test as shown in Fig. 4, Table.S2, and Fig S2. The trend rate for the initial 3 years was excluded both for precipitation and discharge because it is very difficult to examine an upward/downward trend with a small moving window. The value of MK-test significance is estimated on the window of increasing width from 15 years up to the entire period spanned by series and running from the start to the end of the record. The long-term annual precipitation series exhibited varying interdecadal rising and falling trends. For the annual scale, the trend analysis revealed a significant increase in precipitation for Shyok, Shigar, and UIB Kharmong basins and slightly significant for Indus Downstream and Gilgit basins. Although, Hunza and overall flow at Besham Qila also depict an increasing rate of precipitation but is non-significant whereas Astore exhibited a non-significance falling rate in precipitation. The Pettitt test marked 2004 as a change point year, i.e., shift from drying to wetting phase for Indus Downstream, UIB Kharmong-Shyok, and Shingo basins, 2005, was noted for UIB Kharmong. In case of Hunza Gilgit and Shyok basins, a change point has been observed during 2009 and 2012, respectively. In nutshell, precipitation is increasing but nonsignificantly and nonuniformly. It is also noted that Hunza and Shyok basins only exhibited significant change point year by Pettitt test, while all other basins depicted nonsignificant change point year. UIB is characterized by various climatic regimes like westerly disturbances and monsoonal effect orographic disturbance from Tibetan Plateau that makes climatology of this region complex, nonuniform distribution, and results in inconsistency in precipitation anomaly .
On a seasonal scale, trend analysis indicated a significant increasing trend for Hunza, Shyok, UIB Kharmong, and UIB Kharmong-Shyok and a slightly significant rising trend of precipitation for (Gilgit, Shigar) during the winter season. Astore and Indus at Besham Qila expressed a non-significant increasing trend until 2008-2009 followed by a declining rate of precipitation. The 2003 years marked as a change point (drying to wetting phase) for Astore, Hunza (significant), Shingo, and UIB Kharmong basins. Years 2009 and 2011 are noted as 5%-significant change points for Gilgit and Shyok basins. Indus Besham Qila exhibited a change point for the 2013 year but followed inverse phenomena, i.e., wetting to drying phase. For the spring season, UIB Kharmong, Fig. 3 Annual and seasonal precipitation climatology based on CHIRPS gridded datasets corrected using observed data Astore, Shyok, Shigar, and Indus Besham Qila exhibited a non-significant upward trend while, Shingo, Hunza, and UIB Kharmong depicted a non-significant decreasing trend. The change point line was mainly found around the year 2003 for Gilgit, Hunza, Shingo, and UIB Kharmong-Shyok and Shyok basins, 2009, for Shigar and Shyok basins. The increasing rate of precipitation during winter and falling rate during spring is possibly aligned with alteration in westerly precipitation regimes under climate change. The results of winter wetting are also consistent with some previous studies (Cannon et al. 2015;Hasson et al. 2017;Ridley et al. 2013) who also supported rising rate in winter precipitation and drying days during the spring season under climate change is mostly linked with the incursion of westerly precipitation regimes and northward transfer of rainstorm trajectories in UIB.
For the summer season, a non-significant declining trend in precipitation is observed in Shingo/Zanksar, Indus Besham Qila Astore, and UIB Kharmong-Shyok basins. Such fragile monsoon impacts on the lower side of basins are the possible causes of dryness in the summer season. Although precipitation trends are increasing gradually in the remaining basins, they were nonsignificant. Gilgit and Hunza also experienced dryness from 1995 to 2010 that supports the findings (Hasson et al. 2017) regarding precipitation decreased between 1995 and 2012. The change point (drying to wetting phase) varies from 2005 to 2014 for all basins except Shingo that follows in the inverse direction. For the autumn season, Shigar and Gilgit basins revealed a significant upward trend in precipitation. Similarly, Hunza, Shyok, and UIB-Kharmong also show a slightly significant rising trend. Although Indus Besham Qila and Shingo also reveal a rising trend in precipitation, they were nonsignificant. Pettit test found significant change point year for Shyok, Shigar, Gilgit, and Hunza basins while the rest of basins behaved nonsignificantly.

Long-term temperature trends
The trends of long-term annual and seasonal minimum, maximum, average temperature, and diurnal temperature range (DTR) (T min, T max, T avg , and T max -T min ) were observed in order to understand precipitation and streamflow behavior. The MK test and Theil-Sen slope were used to check the significance of trends and their slope. Table 6 shows the annual and seasonal slope values at different stations with bold values indicating the trend's significance. The overall main features of minimum temperature consisted of warming during winter and spring while significant cooling was observed during the summer season. The autumn and annual periods experienced a mixed response. In case of maximum temperature, an overall significant increasing rate of temperature was noted for winter, spring, and annual season while a significant cooling was noted for the summer season. Maximum and minimum winter temperature presented more warming trends than annual time series. Similarly, average temperature also followed significant warming during winter, spring, and annual season while a significant cooling is observed in the summer season. The high agreement of an upward trend for maximum and the average temperature is also associated with the increasing streamflow for all stations during the winter and spring seasons as shown in Fig. 5. DTR also generally displays a significantly increased rate of temperature both for seasonal and annual scales except for the Chilas station.

Variability and trends of the discharge in Upper Indus Basin
The variability and trends of the river discharge were evaluated in different sub-basins as shown in Fig. 5, Table S1, and Fig.S3. The trend magnitude is represented with upward and downward triangles only starting with 20 years time window just for reasons of clarity of the figure. The annual streamflow trends expressed a strong significant increasing trend at Indus Kachura and slightly significant at Shyok and Astore stations. Although streamflow is also rising for Indus (Besham Qila), Bunji, and UIB-Kharmong stations, the rates of increase were not statistically significant. In contrast, streamflow followed non-significant downward trends for Indus at Kharmong, Besham Qila, and Yogo stations. Change point analysis with Pettitt test was also studied for respective streamflow stations. The results marked the year 2004 for Indus at Besham Qila, UIB at Kharmong, and UIB at Kharmong and Shyok, and the year 1988 and 1994 as change point (drying to wetting phase) for Kachura and Shyok stations, respectively. The Pettitt test also depicted a 5%-significant change point in streamflow for Indus (Kachura and Bunji) and Astore and Shyok stations.
In the winter season, most of the subbasins showed a significant increase in streamflow. Although Shyok and Kharmong stations also followed positive trends in streamflow, they were not statistically significant. The increasing streamflow trends in the winter season are consistent with the increase of winter precipitation observed in most catchments, being significant in Shyok, Hunza, and Shigar and also with earlier studies (Khattak et al. 2011;You et al. 2017) that reported climate warming causing early snowmelt. Our analysis of temperature also reported similar results as shown in Table 6. Similarly, the change point analysis also confirms a 5%-significance change at Besham Qila, Kachura, Shyok, Bunji, and Astore in Table S1.
In case of the spring season (MAM) streamflow overall shows an increasing trend in line with the winter season and this can be explained by temperature warming and resulting in earlier and more intense snow-and ice-melt. Indus at Bunji and Astore revealed a statistically significant increase in streamflow. The slightly rising trend in streamflow was also observed in Indus Downstream at Besham Qila, Kachura, Shyok, UIB-Kharmong, and UIB-Kharmong and Yogo. Moreover, change point analysis exhibited significant changes in streamflow for Besham Qila, Astore, and Bunji shown in Fig. 5   For the summer season (JJA), the trend analysis shows a 5%-significance increase in streamflow for Indus at Kachura. The change point analysis also confirms significant variability in streamflow at Kachura. Similarly, slightly significant upward trends were also seen at Astore station. In contrast, a significant decrease in streamflow was observed at Kharmong station. Although Indus at Behsam Qila and Bunji and Shyok also uncovered a declining trend in streamflow, these were not significant. Such a long-term decrease in discharge behaviors is consistent with some previous studies for Indus at Besham Qila and Kharmong by (Arfan et al. 2019;Yaseen et al. 2020) and Indus at Kachura by (Farhan et al. 2015). Similarly, trend behavior for Hunza and Shyok sub-basins is also in agreement with (Mukhopadhyay and Khan 2015). The decrease of flow during the summer season can also be associated with declining temperatures.
These finding demands a serious attention from policy makers and other stakeholder agencies because the variability of summer runoff significantly affects water availability in downstream areas in Indus Basin. In fact, about 70-75% of the Indus flow is generated during the summer season. Such changes in flow trends result in a significant reduction in water availability expected in the coming years. As a major share of this water is being used in the agricultural sector during the summer and winter season in Pakistan, if this trend of flow continues gradually in the long term, the reservoirs, farming, and other water resource management operations must also be implemented and need to adapt accordingly.
In the autumn season (SON), streamflow exhibited a slightly significant rising trend for Indus at Kachura, Shyok, Astore, Bunji, and Besham Qila stations. Similarly, Kharmong station showed a non-significant rising trend in streamflow whereas Besham Qila-Kharmong and Besham Qila-Kharmong-Yogo depicted a slightly significant declining trend in streamflow. The change point analysis marked significant variation in streamflow for Bunji and Astore rivers. The finding of this study is also consistent with a recent study (Yaseen et al. 2020).

Rainfall-runoff relationship in the Upper Indus Basin
The rainfall-runoff relationship and runoff coefficients are the first-order representations of under-or overestimation of precipitation in the watershed. The annual and monthly runoff coefficients at the sub-basin scale were developed as shown in Table 7 and Fig.S4. The results show that both precipitation datasets (CHIRPS and ERA5) including observed values are not able to close the water balance because runoff coefficients (Q/P) higher than one have been calculated in the majority of sub-basins. Higher values of runoff coefficient for Gilgit, Shigar, and Shyok basins depict negative mass balance in these basins. Similarly, Astore, Hunza, and Besham Qila exhibited higher precipitation values than river discharge except in the summer season, because of snow and a glacial melt and have a natural to negative balance at the monthly scale. The Indus Downstream with outlet point at Besham Qila which merges drainage of all upstream sub-basin experiences positive to slight negative mass balance during the summer season. These results need to search for possible explanations.
Overall, ERA5 performed better compared to the CHIRPS precipitation dataset for closing the water balance. Although the CHIRPS dataset has good agreement with observed precipitation values, it is still unable to close the water balance in the majority of basins. Immerzeel et al. (2015) suggested that reanalysis products based on the ECMWF IFS forecast model such as ERA5 can be used to validate atmospheric convergence. ERA5 incorporated fully coupled components of atmosphere, land surface and ocean waves that are useful for closing atmospheric water balance. Gao and Liu (2013)  argued that mountain regions exhibited higher values of runoff coefficient due to greater magnitude of surface runoff, shallow soils, steep slopes, permafrost, and glaciers. In case of UIB, previous studies in this region and neighborhood glacierized catchments also indicated runoff coefficient values greater than one (Adnan et al. 2017;Dahri et al. 2016;Immerzeel et al. 2015;Wortmann et al. 2018). Siddique and Hashmi (2012) found 10%, 25%, and 65% contribution of rainfall, snowmelt, and glacier in the annual flows of Indus River at Tarbela outlet. The results of these studies about higher values of runoff coefficient also support our results. However, values are slightly changed from the current study due to the use of different gridded datasets, size, and location of study area and in any case values higher than one are hardly acceptable.
Generally, it is believed that higher values of runoff coefficient indicate glaciers retreat and alteration of catchment hydrology. However, it is not only the single possible reason for negative mass balance. There are some other factors such as under catch observed precipitation as well as the production mechanism of various gridded precipitation as reported in UIB (Immerzeel et al. 2015;Kääb et al. 2012), and it is also evident in our results as shown in Table 4 and Fig.S3. The discharge values greater than precipitation in Shigar, Shyok, and Gilgit basins might also be associated with under catch precipitation due to the non-availability of observed gauge stations at high elevation in UIB, where the orographic effect on enhancing precipitation could be relevant and because of possible systematic errors in measuring solid precipitation (Eccel et al. 2012). Similarly, some other mass balance studies (Brun et al. 2016;Gao and Liu 2013) also reinforce our conclusion that glacier retreat is only a partial reason for the missing water volumes in UIB for closing the water balance.

Discussion
The diverse climatic signals and contrasting hydrological regimes observed in the UIB are the main sources of the uncertainties affecting the assessment of the key components of the hydrological balance, as precipitation, snow, and ice accumulation and melt and runoff. A clear example of such an inconsistent behavior is the difference between accumulation patterns based on various remote sensing data acquisition techniques and the geodetic mass balance as reported in multiple studies (Immerzeel et al. 2015;Krakauer et al. 2019;Lutz et al. 2016b).
Based on the results of the analysis, it is concluded that the CHIRPS dataset performs well with respect to observed gauge precipitation with the lowest BIAS, MAPE, and RMSE at the annual and seasonal scales, followed by ERA5. The basin-wide corrected monthly precipitation values from CHIRPS and ERA5 and their corresponding runoff coefficient values from each sub-basin are illustrated in Fig.S4. The results of rainfall-runoff, based on a novel combination of gridded datasets and comprehensive ground observations, are in good agreement with some previous studies (Dahri et al. 2016;Kääb et al. 2015). The higher values of runoff in Gilgit, Shigar, and Shyok imply a significant contribution of glacier retreat and snowmelt, as well as undercatch precipitation. It is also concluded that ERA5 precipitation proved to be a better dataset in terms of closure of the water balance.
In the second part of this study, varying positive and negative trends for both precipitation and runoff at seasonal and annual scales in all sub-basins are reported. Previous knowledge about the hydro-climatic trend is mainly confined up to linear trend analysis or with specific time intervals, not explaining non-stationary precipitation and discharge variability within the decadal to interdecadal time scale. The reliable knowledge about hydro-climatic variability over the UIB is very challenging for effective management and precise usage of available water resources in downstream areas (Hasson et al. 2017). In summary, precipitation exhibited greater seasonal than annual variations. Although precipitation is increased annually, its behavior is non-significant except for Shyok and Shigar basins. Pettitt test indicates that change points (drying to wetting phase) mostly lie annually from 2005 to 2010 in the majority of the basins.
An overall increasing trend of winter precipitation is found in all sub-basins. Such a rising rate of precipitation can be due to a significant contribution of winter westerlies regimes and a transfer of rainstorm trajectories in UIB. The results of higher rates in winter precipitation also consensuses with previous studies (Krakauer et al. 2019;Latif et al. 2018;Yaseen et al. 2020). In spring, the majority of glacierized catchments show a downward trend in precipitation. On the other hand, Indus Besham Qila, Astore, and UIB Kharmong indicate the increasing rate of precipitation, but they are statically not significant. Change point analysis also did not record well any transition phase (drying to wetting) in all sub-basins. In summary, spring is drying, as is also reported in some recent studies (Yaseen et al. 2020).
In summer, the basin is not showing any significant trend in the precipitation amount. On the other hand, some basins (Shyok, Gilgit, Hunza, and Shigar) show a rising rate of precipitation, but none of them is statistically significant. The results showing a decreasing rate of summer precipitation align with previous studies (Cannon et al. 2015;Latif et al. 2018;Rizwan et al. 2019). Lutz et al. (2016a) found a clear shift of the summer long-term rising precipitation trends to drying, revealing a transition towards weaker monsoonal influence at lower levels. In order to crosscheck this hypothesis, it would be better to analyze seasonality in precipitation and streamflow by modeling meltwater runoff in the selected area under different climatic conditions. It will be discussed in the future perspective of this study. The Indus Downstream with outlet point at Besham Qila usually receives 70% of the annual precipitation in the summer season. This water is stored in two major reservoirs, Tarbela and Mangla, for the next cropping season, known as Kharif and Rabi season, when rice and wheat are cultivated in major downstream areas of the UIB. If the same downfall trend of precipitation continues in the future, it will reduce water availability, ultimately putting further stress on the already dwindling water reserves of Pakistan.
Concerning streamflow, variabilities are more pronounced seasonally than annually. Results indicate that winter and spring streamflow discharge significantly or slightly significantly increased in all sub-basins, whereas it decreased in summer. Yaseen et al. (2020) suggested that a rising trend of winter discharge is mainly linked with westerly precipitation regimes because a major portion of UIB hydrology is dominated by westerly disturbances rather than monsoon offshoots. There are also different significant interpretations about these flow dynamics. One reason could be found in the significant warming in winter and spring, as shown in Table 6, whereas summer cooling caused early snowmelt during spring and less flow available during summer (i.e., decreasing trends in summer discharge show lower melting rates in summer, resulting in potential stability of glaciers and consequently positive basin storage).

Conclusions
The study presents a comprehensive hydro-climatic trend and precipitation anomaly analysis for the UIB at the subbasin scale. The primary objective of this study is to evaluate the performance of four gridded precipitation datasets for developing precipitation climatology and check its reliability for the UIB. The datasets were examined for an overlapping period spanning from 1995 to 2017 at the seasonal and at the annual scale. Based on the results, it is found that the performance of the CHIRPS dataset is good to describe the distribution of observed precipitation with the lowest BIAS, MAE, RMSE, and MAPE, followed by ERA5. The mean annual corrected precipitation was calculated as 536 mm/year in the UIB gauged at Besham Qila. The precipitation climatology exhibited a higher rate of precipitation in the lower part of the basin for both the annual and the seasonal scales. The runoff coefficient for CHIRPS and ERA5 is though greater than one in some basins, making the water balance unrealistic. There can be two main reasons: (1) underestimate of precipitation, as most of the monitoring stations in the UIB are valley-based and do not represent the true basin hydrology in the high elevation bands and (2) glacier retreat and early snowmelt due to global warming and elevation-dependent warming. However, there are small chances for glacier retreat because most glaciers, especially in the Karakorum, have been advancing or in stable conditions in the last decade (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) (Berthier and Brun 2019). Meanwhile, the precipitation rate declines with elevation annually, rises during the winter and spring seasons but decreases during the summer season. These issues demand further investigation, as they are affecting the contribution of glaciers and snowmelt in total flow from each sub-basin. The findings of this study would be helpful to understand the discrepancies between the observed and the gridded precipitation datasets referring to the UIB and may have a substantial impact on studies related to the designing, planning, modeling, and management of the water resources under climate change. The results of the study also recommend that gridded precipitation should be corrected before its usage in hydrological modeling studies, especially in those involving glacierized catchments. The anomaly method proved to be worthwhile for assessing precipitation climatology, especially in data-scarce regions with a sparse monitoring network.
In the second part of this study, annual and seasonal precipitation revealed significant variability seasonally rather than annually. Summer is drying, while winter is wetting. The increasing rate of precipitation was also seen during spring in some basins, but they were not statistically significant. Similarly, trend analysis of observed streamflow at various gauge stations in the UIB facilitates understanding of comprehensive water balance for the region. Like precipitation variability, the streamflow one is more pronounced seasonally rather than annually. At the annual scale, trend analysis of discharge shows a slightly significant increasing trend at the Indus River Kachura, Shyok, and Gilgat stations, while nonsignificant decreasing trends at Kharmong stations, Besham Qila-Kharmong, and Besham Qila-Kharmong-Shyok stations are found. Seasonal flow analysis reveals a more complex regime: winter (December-February) and spring (March-May) exhibit a rising trend in streamflow, while summer (June-August) shows a declining trend. The seasonal analysis also shows an increasing rate of warming in spring and early seasonal melt discharge from most of the sub-basins, whereas field significant low flow/drying was observed during summer.
The findings of this hydro-climatic analysis are expected to support future sustainable development projects in the study area. For instance, it would be helpful to assist engineers, the government and its organizations, as well as other stakeholder agencies, to set up structural and non-structural measures to handle extreme flood and other natural hazard events, such as building dams and other control structures, lining canals, and watercourse and adopt precision agricultural techniques (drip and sprinkler irrigation). It would also be viable to bridge the gap in terms of water availability and supplies especially in the lower area of the basin, where a major share of this water is being used for growing crops.
These results would facilitate farmers and other stakeholder agencies to set cropping patterns according to water availability under prevailing climatic conditions.