Estimation of the River Flow Synchronicity in the Upper Indus River Basin Using Copula Functions

: In this study, on the basis of the maximum and mean annual values of ﬂows, dependencies between ﬂows recorded in seven water gauges located in the upper part of the Indus River Basin (IRB) in Pakistan were analyzed. First, the non-parametric Mann–Kendall (M–K) test was used to detect trends in the ﬂows. Next, the Pearson’s correlation coe ﬃ cient was applied. Then, the selected copulas were used to ﬁnd joint distributions of the studied time series. In the next stage, the degrees of synchronous and asynchronous occurrence of, respectively, the annual maximum (AMAXF) and mean annual ﬂows (MAF) were calculated. The study revealed that correlations between the ﬂows in selected gauge stations were very strong and statistically signiﬁcant. These results were conﬁrmed by the synchronicity analysis carried out with the help of the copula functions. The highest relationship was detected in the case of gauges Besham Qila and Kachura on the Indus mainstream, while the lowest was detected in gauges Besham Qila and Naltar on the Naltar River. These ﬁndings can be of high practical value in the ﬁeld of sustainable water resource management, including for ﬂood protection, agricultural water supply, reservoir water storage, and hydropower generation in the IRB.


Introduction
Sustainable management of water resources is of particular importance in water-stressed and populous areas. Pakistan, with nearly 213 million people, ranks fifths in terms of the world's population. Projections suggest that Pakistan's population will exceed 300 million by 2047, driving water demands much higher. Pakistan's water resources are already highly stressed and will become increasingly so with projected population changes [1]. Without serious demand management and reform, and if the climate warms rapidly, water demand could increase by nearly 60 percent by 2047 [2].
Pakistan's economy heavily relies on the agricultural sector, which accounts for a quarter of its GDP and employs two-fifths of total labor force. Pakistan is located in a semi-arid to arid region, where rainfall is very low. The agriculture mainly depends on the Indus River System (IRS) for 90% of its irrigation needs, which also supplies by about 30% of its total energy generation [3].
The Indus River is the largest river of Pakistan. Its total length is 3180 km, and basin area 1,165,000 km 2 . It originates in the Tibetan Plateau, runs south, and discharges into the Arabian Sea. Its annual runoff is around 207 km 3 , which makes the Indus one of the largest rivers in the world in terms of annual flow. The Indus River Basin (IRB) covers 65% of Pakistan and represents over 95% of its water resources. The Indus River and its tributaries provide the main source of water for Pakistan's vast irrigation network and the hydropower generation for its power grid in an otherwise that the vast majority of the water systems in the UIB were seriously impacted by environmental change. Yassen et al. [24] analyzed hydrological trends and their variability in the Mangla River watershed in 1971-2010, and revealed that the maximum flows for the spring, summer, and autumn seasons had decreasing trends in all sub-basins. Hasson et al. [25] drew a conclusion that the hydrology of the UIB was dominated with the meltwater runoff, which ensured the crucial water supply to the largest reservoir in Pakistan for reducing the ongoing electric shortfall by its use for hydropower generation, and contributing to the economy through its use mostly for irrigated agricultural production downstream. Shamshad et al. [26] analyzed the impact of climate change on catchments of the Indus, Hunza, and Jhelum River basins, and determined an increasing trend of winter warming and summer cooling in the studied basins. Faiz et al. [27] applied the Shannon information entropy theory to assess precipitation variability and uncertainty of stream flow in the Hindu Kush, Himalayan, and Karakoram River basins of Pakistan, and found that the rivers with high stream flow variability also showed low entropy of their distribution and, therefore, higher stream flow concentration in the annual cycle. The study of [28] on the spatial and temporal dynamics of precipitation in agricultural zones of arid areas of south-west Pakistan indicated that climate change was going to seriously affect the region, as a decreasing trend prevailed in most of the cases and stations.
As can be seen from the above review, different analyses of the hydrology and water resource management of the UIB have been carried out, yet the synchronous and asynchronous manifestations of the relations of both the annual maximum flow (AMAXF) and mean annual flow (MAF) in these research areas have not been examined until now. While being aware of the importance of the Indus River's tributaries in the formation of water resources used for the socio-economic purposes, it seems equally important to recognize relationships between, respectively, the AMAXF and MAF of the Indus River and its tributaries in terms of their synchronous occurrence. In order to overcome that research gap, this study was designed for better understanding of the synchronicity and asynchronicity of the river flows of the river systems of the Himalaya and Hindu Kush-Karakoram Mountains in Northern Pakistan. This was achieved with the help of the copula function. This paper is the first attempt to determine the respective dependencies between AMAXF and MAF in the UIB in Pakistan using selected Archimedean copulas.
Studying runoff variation is the assurance for water resources' integrated development, scientific management, and optimal operation [29]. The synchronicity and asynchronicity of the AMAXF and MAF and the probability of their occurrence provide important information about the hydrologic activity of the basin, including transformations of the flood wave and their overlapping in the IRB. Thus, their recognition is of key importance for the economy of Pakistan. According to [30], the Indus River Basin is a very flood-prone area. Rising population pressures, climate change, and a continuous degradation of ecosystem services have resulted in increased flood risks, which are further exacerbated by inadequate flood planning and management. A total of 21 major floods occurred between 1950 and 2010 in the Indus Basin, causing cumulative direct economic losses of about 19 billion USD (in 2010 dollars), killing 8887 people, and damaging or destroying a total of 109,822 villages (within an area of around 446,000 km 2 ). The catastrophic floods of 2010 in the Indus Basin were the largest in recorded history [31], and affected more than 14 million people in Pakistan [32].
The aim of this study is to analyze the probabilities of, respectively, synchronous and asynchronous occurrence of the AMAXF and MAF of rivers in the UIB in Northern Pakistan. Apart from the purely scientific and theoretical outcomes of the study, the obtained results can be of particular importance for river basin management, flood protection, and sustainable utilization of its water resources. Recognition of such dependencies can also be helpful in proper operation of the Tarbela Dam, which collects waters of the analyzed rivers, or the Diamer-Bhasha and Mohmand dams in the UIB, which are currently being constructed.

Study Area
The UIB is a transboundary basin, shared by four countries: Afghanistan, China (Tibet), India, and Pakistan. Its total area in the gauge Besham Qila in Pakistan is estimated at around 164,800 km 2 [8]. Precipitation is not uniform over the UIB. Annual precipitation ranges between 100 mm in the Gilgit area to 1000 mm in Besham Qila, and a maximum of 1500 mm on the mountain slopes at Murree. The highest precipitation is recorded in March and April, while the lowest is in December and January [33]. The hydrological zone of the UIB is situated in the high altitudes of the Karakoram, Himalaya, and Hindu Kush ranges. Snowfall at higher altitudes accounts for most of the river runoff. Approximately three-fourths of the UIB's discharge are estimated to be a result of melting snow and ice, while summer rainfall plays a smaller role in contributing to the total discharge [13].
The study is based on the values of AMAXF and MAF collected in seven gauges situated in the UIB in Pakistan: Kachura and Besham Qila on the Indus River mainstream, Yogo on the Shyok River, Naltar on the Naltar River, Dainyor on the Hunza River, Alam Bridge on the Gilgit River, and Doyian on the Astore River ( Figure 1) in 1974-2004.

Study Area
The UIB is a transboundary basin, shared by four countries: Afghanistan, China (Tibet), India, and Pakistan. Its total area in the gauge Besham Qila in Pakistan is estimated at around 164,800 km 2 [8]. Precipitation is not uniform over the UIB. Annual precipitation ranges between 100 mm in the Gilgit area to 1000 mm in Besham Qila, and a maximum of 1500 mm on the mountain slopes at Murree. The highest precipitation is recorded in March and April, while the lowest is in December and January [33]. The hydrological zone of the UIB is situated in the high altitudes of the Karakoram, Himalaya, and Hindu Kush ranges. Snowfall at higher altitudes accounts for most of the river runoff. Approximately three-fourths of the UIB's discharge are estimated to be a result of melting snow and ice, while summer rainfall plays a smaller role in contributing to the total discharge [13].
The study is based on the values of AMAXF and MAF collected in seven gauges situated in the UIB in Pakistan: Kachura and Besham Qila on the Indus River mainstream, Yogo on the Shyok River, Naltar on the Naltar River, Dainyor on the Hunza River, Alam Bridge on the Gilgit River, and Doyian on the Astore River (   Table 1.  The data were obtained from the Surface Water Hydrology Project of the Water and Power Development Authority (WAPDA) of Pakistan. Major characteristics of the analyzed rivers and flows are given in Table 1.
The Astore River lies in the north-facing slopes of the Himalayan Range that are covered with glaciers, near the Burzil Pass. It is fed by small snow-covered streams and joins the Hunza River. The Hunza River discharges into the Gilgit River at Dainyor, and the Gilgit River joins the Indus River near Bunji, while, from the south-east, the Astore River falls into the Indus River just below Bunji [34]. According to the available flow records, the annual river flows gauged at Hunza, Astore, and Gilgit are equivalent to 742 mm, 1084 mm, and 800 mm of water depth, respectively [17].

Mann-Kendall Test
In the first stage of the study, tendencies of fluctuations of the AMAXF and MAF in the seven respective gauges were analyzed. This involved the application of the rank-based non-parametric Mann-Kendall (M-K) test, which detects trend in temporal sequences [35]. The M-K test is applicable in cases when the data values x i of a time series can be assumed to obey the model: where f(t) is a continuous monotonic increasing or decreasing function of time, and the residuals ε i can be assumed to be from the same distribution with zero mean. It is therefore assumed that the variance of the distribution is constant in time.
The M-K test statistic S is calculated using the formula: where x j and x k are the time series observations in chronological order in years j and k, j > k, respectively, n is the length of time series, and: An upward (increasing) or downward (decreasing) trend is determined by a positive or negative value of Z. First, the variance of S is computed by the following Equation (4), which takes into account that ties may be present: where q is the number of tied groups and t p is the number of data values in the pth group. The values of S and VAR(S) are used to compute the test statistic Z as follows: Then, the null hypothesis of no trend, H 0 , is tested in order to accept or reject it. The x i observations are randomly ordered chronologically, contrary to the alternative hypothesis H 1 , where there is an increasing or decreasing monotonic trend. In this study, the null hypothesis (H 0 ) was that there had been no trend in AMAXF and MAF over time. The alternate hypothesis (H 1 ) was that there had been a trend (increasing or decreasing) over time.
The test statistic Z (normal approximation) is computed if all time series are longer than ten. The statistic Z has a normal distribution. The absolute value of Z can be compared to the standard normal cumulative distribution in order to identify if there is a monotone trend or not at the specified level of significance.
A huge number of studies used the M-K non-parametric test because of its ease of usage and flexibility with missing values, but also because of appealing features like skewed distribution, as mentioned in papers of [36,37].
In this study, the M-K test employed the MAKESENS stencil in the Microsoft Excel software, developed by scientists from the Finnish Meteorological Institute [38].

Correlation
Next, the sample Pearson's correlation coefficient was applied in the analysis of correlations of the AMAXF and MAF. The Pearson's correlation coefficient is a measure of the strength of the association between two variables. Then, statistical significance was analyzed at a level of p = 0.05. Given paired data (x i , y i ), . . . , (x n , y n ) consisting of n pairs, r xy is defined as: where n is sample size, x i , y i are the individual sample points indexed with i, and x, y are the sample means.

Application of the Copula Theory
Sklar [39] introduced the concept of the copula and defined it as a joint distribution function of standard uniform random variables. Copulas offer a flexible way of describing the dependence structure of multivariate data independently of their marginal probability distributions. Additionally, copulas are powerful tools for modeling and sampling multivariate data that are nonlinearly interrelated [40]. Numerous successful applications of copula modeling have been made, most notably in survival analysis, actuarial science, and finance [41]. It is found that copula-based flood frequency analysis performs better than conventional flood frequency analysis, as joint distribution based on copula fits the empirical joint distribution (i.e., from observed data using plotting position formula) better than that established by standard joint parametric distribution. De Michele and Salvadori [42] examined that copulas within hydrology are utilized to analyze the severity of the events by studying the flooding as a collaborative behavior of the volume and peak. Bárdossy and Pegram [43] studied that copulas have been applied to depict the spatiotemporal uncertainties of precipitation or the homogeneity parameters of land water [44]. Zhao et al. [45], using the copulas, analyzed the frequency analysis of extreme precipitation in the Pearl River Basin in South China. Kuchment and Demidov [46] applied the copula theory to study the spring flood flows and flow volumes in the Belaya and Vyatka Rivers in Russia. Rizwan et al. [47] used the Archimedean family of copulas to construct the bivariate distribution of flood peaks and volumes of the Indus and Jhelum Rivers in Pakistan. Sugimoto et al. [48] studied hydrological time series using copulas to detect catchment characteristics and anthropogenic impacts in the upper Rhine region in Germany. Zhou et al. [49] used the copulas to evaluate the synchronous and asynchronous runoff-sediment encounter probability in the Dongting Lake Basin, and You et al. [50] carried out the probabilistic analysis of the Weihe and Jinghe Rivers' runoff-sediment characteristics. Bacova, Mitková, and Halmova [51] applied the copulas for modeling of flood peak flows, volume, and duration of the Danube River in Bratislava, Slovakia. Based on selected Archimedean copulas, correlations between the maximum [52] and mean [53] annual water levels in the Polish coastal lakes with the Baltic Sea water levels were assessed. Other examples of copula applications in hydrology and water resource management are given in [54].
In this study, at the beginning, the best-matching statistical distributions were selected for the analyzed data series. The following distributions reflect the maximum and mean values in series of hydrological data: Weibull, Gamma, Gumbel, and log-normal. Parameters of the distributions were estimated by means of the highest probability method. For the purpose of assessment of the goodness of fit of a given distribution to the data series, the Akaike Information Criterion (AIC) was applied [55], calculated from the following formulas: where MSE is the mean square error and N is the sample size, or: The best model is the one that has the minimum AIC value.
In the next step, the copula method was used to construct the joint distribution of flow series in two compared gauge stations. In general, a bivariate Archimedean copula can be defined as: where u and v are marginal distributions, and the θ subscript of copula C is the hidden parameter of the generating function φ. φ is a continuous function, called the generator, and is strictly decreasing In hydrological studies, the Archimedean copula family is often applied. A large variety of copulas belong to this copula family, and they can be used when the correlations among hydrologic variables are positive or negative. The proofs of these properties have been reported in [56,57]. For this reason, the one-parameter Archimedean copulas, including the Clayton family, the Gumbel-Hougaard family, and the Frank family, were used in this study. The Gumbel-Hougaard and Clayton copula families are appropriate only for the positively correlated variables, while the Frank family is appropriate for both negatively and positively correlated variables (Table 2). Table 2. Copula function, parameter space, generating function φ(t), and functional relationship of Kendall's τ θ with copula parameters for selected single-parameter bivariate Archimedean copulas.

Copula Family
where D k (x) is the Debye function; for any positive integer k, The best-fitted joint distribution was selected by comparison with the empirical joint distribution using the Akaike Information Criterion (AIC), as mentioned earlier. For each compared pair of series, based on previously calculated parameters of statistical distribution, 5000 hypothetical points were randomly generated. They were used for the selection of the best fitted family of copulas for a given pair of series, and then for the development of an appropriate copula. Based on empirical pairs of values for particular years and generated hypothetical points, graphs with probability curves (expressed in return periods) were developed- Figures 2-4. The next stage involved the calculation of the degree of synchronicity (synchronous occurrence) and asynchronicity (asynchronous occurrence) of AMAXF and MAF in the UIB. For each pair of gauges, probability curves at the level of 62.5% (return period once in 1.6 years), 37.5% (once in approximately 2.7 years), 20% (once in 5 years), 10% (once in 10 years), 2% (once in 50 years), 1% (once in 100 years), 0.5% (once in 200 years), and 0.2% (once in 500 years) were presented, respectively ( Figure 2).
The obtained data were then analyzed based on probabilities of 62.5% and 37.5% [58,59]. Nine sectors were designated, representing different relations between probable AMAXF and MAF. On the basis of generated points with a distribution imitating shared distribution of values from compared water gauge stations and their participation in particular sectors (Figure 2), three sectors with synchronous occurrence of AMAXF and MAF were designated: The next stage involved the calculation of the degree of synchronicity (synchronous occurrence) and asynchronicity (asynchronous occurrence) of AMAXF and MAF in the UIB. For each pair of gauges, probability curves at the level of 62.5% (return period once in 1.6 years), 37.5% (once in approximately 2.7 years), 20% (once in 5 years), 10% (once in 10 years), 2% (once in 50 years), 1% (once in 100 years), 0.5% (once in 200 years), and 0.2% (once in 500 years) were presented, respectively ( Figure 2).
The obtained data were then analyzed based on probabilities of 62.5% and 37.5% [58,59]. Nine sectors were designated, representing different relations between probable AMAXF and MAF. On the basis of generated points with a distribution imitating shared distribution of values from compared water gauge stations and their participation in particular sectors (Figure 2), three sectors with synchronous occurrence of AMAXF and MAF were designated: where: X-values of x-coordinates of generated points; Y-values of y-coordinates of generated points; A 62.5% -value of AMAXF/MAF with probability of exceedance of 62.5%; A 37.5% -value of AMAXF/MAF with probability of exceedance of 37.5%; B 62.5% -value of AMAXF/MAF with probability of exceedance of 62.5%; B 37.5% -value of AMAXF/MAF with probability of exceedance of 37.5%; F-flow; H-high (maximum); M-mean.
The percent contribution of points included in sectors 1, 5, and 9 allowed determination of the degree of synchronicity of AMAXF and MAF between the two analyzed gauges in a given time unit.
Synchronous and asynchronous occurrences of AMAXF and MAF were calculated through the determination of threshold values of probability ranges: • Probable AMAXF and MAF with probability of occurrence of <62.5% were designated as LHF/LMF, • Probable AMAXF and MAF with probability of occurrence in a range of >62.5% and <37.5% were designated as MHF/MMF, • Probable AMAXF and MAF with probability of occurrence of >37.5% were designated as HHF/HMF.
For example, the occurrence of LHF in a given gauge is a synchronous event if LHF also occurs in the other gauge in a given time unit, and asynchronous if MHF or HHF is recorded there.
The sum of the percent shares of synchronous and asynchronous events is always 100%. The mathematical and statistical processing of analysis results employed statistical procedures included in the following software: Excel (Microsoft, Redmond, USA), Statistica (TIBCO Software Inc., Palo Alto, USA), and RStudio (Boston, USA). The implementation of the graphic form employed QGIS (3.6.2. Noosa) and Publisher (Microsoft, Redmond, USA) software.

Mann-Kendall Test
Multi-annual fluctuations of the AMAXF and MAF showed statistically significant increasing trends in two gauges: Naltar on the Naltar River and Kachura on the Indus River. In the rest of the analyzed gauges, a mix of non-significant increasing and decreasing trends was detected. AMAXF in gauge Doyian on the Astore River had increasing trends, while MAF was decreasing. In gauges Yogo on the Shyok River and Alam Bridge on the Gilgit River, AMAXF and MAF revealed decreasing trends, while in gauge Dainyor on the Hunza River, AMAXF had a decreasing trend, and the mean was increasing. In gauge Besham Qila on the Indus River, both AMAXF and MAF showed non-significant increasing trends (Table 3). Statistically significant trends at *** p < 0.001, ** p < 0.01, and * p < 0.05.

Annual Maximum Flows (AMAXF)
The correlation coefficients and the degree of synchronicity of occurrence of AMAXF in the analyzed gauges were determined. The obtained results show different degrees of positive correlations between AMAXF in the analyzed gauges ( Table 4). The strongest (0.66) and statistically significant (p < 0.05) correlation was detected for the pair Besham Qila-Kachura, both located on the Indus River mainstream. Statistically significant were also correlations between the pairs Besham Qila-Yogo on the Shyok River (0.52) and Besham Qila-Doyian on the Astore River (0.48). For the rest of the gauges (Naltar, Dainyor, and Alam Bridge) the correlations with gauge Besham Qila were statistically non-significant and relatively weak. They were 0.31 for gauge Alam Bridge on the Gilgit River and 0.25 for gauge Dainyor on the Hunza River, respectively. The weakest correlation (0.24) was determined for gauge Naltar on the Naltar River.

Mean Annual Flows (MAF)
Coefficients of correlation between all pairs of the mean annual flows in the UIB are high and statistically significant (p < 0.05)- Table 5. The strongest (0.76) correlations were detected for the pair Besham Qila-Kachura on the Indus River mainstream and Besham Qila-Gilgit River. On the other hand, the weakest correlations of the MAF were observed for the pair Besham Qila-Naltar (0.44). Table 5. Results of analysis of correlation and synchronicity between the analyzed gauge stations for mean annual flow (MAF).

Annual Maximum Flows (AMAXF)
Like the above described results of correlation coefficients, modeling correlations with the help of the copula function allows us to determine the strength of correlations between flows in analyzed gauges. Moreover, the copula function gives the possibility of quantifiable determination of the probability of occurrence of each analyzed variable separately and the probability of co-occurrence of AMAXF in gauges, as well as the degree of their synchronicity.
In the case of half of the data pairs, the best fit was the Gumbel copula family, and for the second half of the data pairs, it was the Clayton copula family (Figure 3). The Frank copula family was not chosen for any of data pairs as the best fit in synchronicity analysis of AMAXF. Based on Figure 3, presented below, it can be seen, for example, that the cloud of randomly generated points for the pair Besham Qila-Kachura on the Indus River mainstream (where the highest synchronicity was identified) is more concentrated in the figure axis than in the case of the remaining pairs of gauges.
Sustainability 2020, 12, x FOR PEER REVIEW 11 of 18 In the case of half of the data pairs, the best fit was the Gumbel copula family, and for the second half of the data pairs, it was the Clayton copula family (Figure 3). The Frank copula family was not chosen for any of data pairs as the best fit in synchronicity analysis of AMAXF. Based on Figure 3, presented below, it can be seen, for example, that the cloud of randomly generated points for the pair Besham Qila-Kachura on the Indus River mainstream (where the highest synchronicity was identified) is more concentrated in the figure axis than in the case of the remaining pairs of gauges.  The synchronicity of AMAXF in the analyzed gauges is relatively low and oscillates between 42.06% and 61.26%, and asynchronicity, in consequence, between 38.74% and 57.94% (Table 5). Relations between flows resulting from the synchronicity analysis are proved by results obtained in the correlation analysis (Table 5)-the higher value of correlation coefficient, the higher synchronicity level.
The highest synchronicity of AMAXF is observed in gauges Besham Qila and Kachura, both located on the Indus River mainstream (61.26%). This means that AMAXF in gauge Kachura in a hydrological year will be in the same range of probability of occurrence as AMAXF in gauge Besham Qila (i.e., one of three possible synchronous situations will occur: LHF A -LHF B , MHF A -MHF B , HHF A -HHF B ) with a probability of approximately 61%, which statistically means three times in every five years). Given that these gauges are located on the same river at a distance of about 440 km, this is relatively weak level of synchronicity. On this basis it can be concluded that tributaries of the Indus between those two gauges significantly moderate AMAXF of that river.
The lowest synchronicity (and, consequently, the highest asynchronicity) is characteristic for gauges Besham Qila and Naltar (42.06%). Only gauge Yogo on the Shyok River has synchronicity of AMAXF slightly higher than 50%. In case of the other gauges located on tributaries of the Indus River asynchronicities are higher than synchronicities. It means that the flows of the Indus tributaries on the analyzed section are asynchronous in relation to gauge Besham Qila, so they are in a different range of probability in the same hydrological year. Consequently, the probability of occurrence of similar AMAXF in a given hydrological year in the analyzed gauges is lower than the probability of occurrence of asynchronous events. It can mean that no single tributary shapes the flow of the Indus River in Besham Qila, but together they moderate it in significant way, as mentioned earlier.
More detailed analysis of results of the percent share of synchronous and asynchronous occurrences of AMAXF in respective sectors (Table 6) shows that, for example, in three cases (gauges Yogo, Doyian and Kachura) the probability of occurrence of HHF A -HHF B (sector 9, see Figure 2) is higher than LHF A -LHF B (sector 1, see Figure 2). This means that in relation to gauge Besham Qila in these three gauges the synchronous occurrence of extreme floods is slightly more probable than the small ones. Table 6. Percent share of synchronous and asynchronous occurrences of annual maximum flow (AMAXF) in analyzed gauges in nine sectors.

Sector
Besham Each copula family (Gumbel, Frank and Clayton) was best fitted twice in terms of MAF ( Figure 4). Also in this case it can be seen that clouds of generated points are more concentrated in the figure axis in the case of gauge pairs with the highest synchronicity ( Figure 4). The synchronicity of MAF was calculated for all pairs of the investigated gauges- Table 7. It ranges from 49.00% to 66.96% (Table 7). So, synchronous occurrence of the mean annual flows is more probable than AMAXF.
The mean annual flows are asynchronous only in two gauges (Dainyor and Naltar) in relation to Besham Qila (Table 7). The most synchronous (almost 67%) occurrence of MAF characterizes gauge Alam Bridge. It means that about two times every three years, MAFs in Alam Bridge on the Gilgit River are in the same probability range as in Besham Qila on the Indus River. It is even higher than in gauge Kachura (61.82%), and likewise Besham Qila located on the Indus mainstream. This can mean that the Gilgit River has a significant impact on shaping flows on the Indus mainstream. The synchronicity of MAF was calculated for all pairs of the investigated gauges- Table 7. It ranges from 49.00% to 66.96% (Table 7). So, synchronous occurrence of the mean annual flows is more probable than AMAXF.
The mean annual flows are asynchronous only in two gauges (Dainyor and Naltar) in relation to Besham Qila (Table 7). The most synchronous (almost 67%) occurrence of MAF characterizes gauge Alam Bridge. It means that about two times every three years, MAFs in Alam Bridge on the Gilgit River are in the same probability range as in Besham Qila on the Indus River. It is even higher than in gauge Kachura (61.82%), and likewise Besham Qila located on the Indus mainstream. This can mean that the Gilgit River has a significant impact on shaping flows on the Indus mainstream. Further analysis of respective sectors shows that in reference to Besham Qila, for example, in gauges Yogo and Kachura, the probability of synchronous occurrence of HMF A -HMF B (sector 9, see Figure 2) is higher than LHF A -LHF B (sector 1, see Figure 2). This means that in wetter years (when the mean flows are higher than on average), the synchronous occurrence of MAF in these two gauges and in Besham Qila is slightly more probable than in drier years.

Discussion and Conclusions
This research aimed at detection of synchronous and asynchronous occurrence of AMAXF and MAF in the UIB in Northern Pakistan in 1974-2004. The M-K test applied in the first part of the study revealed statistically significant increasing trends for both AMAXF and MAF in two gauges: Naltar on the Naltar River and Kachura on the Indus River mainstream. However, increasing non-significant trends were also detected in gauge Besham Qila on the Indus River mainstream. In the rest of the analyzed gauges, only statistically non-significant decreasing (Yogo on the Shyok River, Alam Bridge on the Gilgit River) or alternately decreasing and increasing (Dainyor on the Hunza River, Doyian on the Astore River) trends for AMAXF and MAF were found.
In the second part of the study, positive correlations, yet of different degrees, between AMAXF and MAF in the respective gauges were determined. While for AMAXF, the highest (0.66) statistically significant (p < 0.05) correlation was found for the pair Besham Qila-Kachura, both on the Indus River mainstream, for MAF, the highest (0.76) statistically significant correlation was detected for two pairs of gauges, namely Besham Qila-Kachura and Besham Qila-Alam Bridge on the Gilgit River. In the case of AMAXF, relatively high and statistically significant (p < 0.05) correlations were detected for the pairs Besham Qila-Yogo on the Shyok River (0.52) and Besham Qila-Doyian (0.48) on the Astore River. In turn, the relationships between gauges Besham Qila on the Indus mainstream and, respectively, Alam Bridge on the Gilgit River (0.31), Dainyor on the Hunza River (0.25), and Naltar on the Naltar River (0.24) were the weakest and statistically non-significant. Contrary to them, coefficients of correlations calculated for MAF were statistically significant for all pairs of gauges, with the weakest relationships for Besham Qila-Dainyor on Hunza River (0.47) and Besham Qila-Naltar on the Naltar River (0.44). These results proved higher correlations of MAF than AMAXF.
The copula function allowed determination of the synchronous and asynchronous encounter probability and the degree of synchronicity of, respectively, AMAXF and MAF in the respective gauges in relation to the closing profile of Besham Qila. It was found that the degree of synchronicity corresponds to the strength of relationships measured by the coefficient of correlation. The synchronicity of AMAXF is relatively lower compared to that of MAF. For AMAXF, the highest synchronicity (61.26%) was found for the pair Besham Qila-Kachura on the Indus mainstream, while for MAF, that pair ranked second (61.82%) after Besham Qila-Alam Bridge on the Gilgit River (66.96%). In turn, the lowest degree of synchronicity for AMAXF and MAF was detected for the pair Besham Qila-Naltar on the Naltar River-42.06% and 49.00%, respectively. These mean that flows on the Gilgit River have the highest relative impact on shaping flows on the Indus River, while these play a smaller role on the Naltar River, which is justified given their flow volume (Table 1). Generally, it can be concluded that while no single tributary shapes flows of the Indus River in Besham Qila, altogether, they moderate them significantly.
More detailed analysis of the obtained results of the copula modeling allowed determination of the percent shares of synchronous and asynchronous occurrences of AMAXF and MAF. On this basis, ranges of probabilities of the synchronous occurrence of extreme floods in respective gauges were calculated. This allowed the conclusion that, for example, with regard to AMAXF, in relation to gauge Besham Qila, in gauges Yogo on the Shyok River, Doyian on the Astore River, and Kachura on the Indus River, the probability of synchronous occurrence of extreme floods is higher than moderate floods, while regarding MAF, their synchronous occurrence in gauges Yogo and Kachura in relation to Besham Qila in wetter years is more probable than in drier years.
These interesting results allow an analysis of the spatial distribution of the synchronous and asynchronous occurrences of AMAXF and MAF in the respective tributaries along the Indus River mainstream. In the case of AMAXF, a relatively high (50.66%) synchronicity in relation to gauge Besham Qila on the Indus mainstream is recorded in gauge Yogo on the Shyok River, but it is significantly lower in gauge Naltar on the Naltar River (42.06%), a tributary of the Hunza River (43.22% in gauge Dainyor), which then discharges into the Gilgit River (42.86% in gauge Alam Bridge). The synchronicity is relatively higher (49.04%) in gauge Doyian on the Astore River. In the case of MAF, the synchronicity is also relatively high (52.62%) in gauge Yogo on the Shyok River, while it is lower in the sub-basins of the Naltar (49.00% in gauge Naltar) and Hunza (49.06% in gauge Dainyor) rivers. However, in gauge Alam Bridge on the Gilgit River, the synchronicity of MAF reaches the relatively highest (66.96%) degree among all the analyzed gauges, while it reaches 59.70% in gauge Doyian on the Astore River. Generally, it should be noted that the Naltar River, which has the lowest rank in the bottom-up hydrographical ordering among the analyzed rivers, also has the highest asynchronicity in relation to gauge Besham Qila on the Indus mainstream. It is also worth noting that the Gilgit River, the first-order tributary of the Indus River, has the highest degree of synchronicity in terms of MAF; however, it behaves asynchronously if one considers AMAXF. In other words, despite the fact that the MAFs on the Gilgit River and the Indus mainstream often fall within the same range of probability, they look different in the case of AMAXF on these two rivers. This can mean that flood waves on the Gilgit River appear regardless of the flood wave on the Indus River mainstream. It can be assumed that this spatial differentiation of the synchronous and asynchronous occurrence of AMAXF and MAF in the respective sub-basins is connected with their hierarchical position within the river system, geographical location, morphology, rainfall, snow cover, and other factors. This needs further and more detailed studies.
Pakistan faces a tremendous shortage of water due to the growing population, increasing agricultural water demand, and decline in reservoir water storage capacity. On the other hand, catastrophic floods frequently affect the Indus River basin, including its upper part. On average, the annual flood damage from 1960 to 2011 was about 1% of the mean annual GDP [30].
The obtained results not only broaden our knowledge about the hydrology of the UIB, but also have a considerable applicative potential. Forecasting AMAXF and MAF of the Indus River with the use of the copula function can help to formulate adaptation strategies and reconsider existing regional development plans, including management projects for irrigation schemes, flood protection, hydropower generation, and sustainable water resources planning in that area and in these regions of Pakistan, which benefit from the UIB's waters.