Non-Stationary Bayesian Modeling of Annual Maximum Floods in a Changing Environment and Implications for Flood Management in the Kabul River Basin, Pakistan

: Recent evidence of regional climate change associated with the intensiﬁcation of human activities has led hydrologists to study a ﬂood regime in a non-stationarity context. This study utilized a Bayesian framework with informed priors on shape parameter for a generalized extreme value (GEV) model for the estimation of design ﬂood quantiles for “at site analysis” in a changing environment, and discussed its implications for ﬂood management in the Kabul River basin (KRB), Pakistan. Initially, 29 study sites in the KRB were used to evaluate the annual maximum ﬂood regime by applying the Mann–Kendall test. Stationary (without trend) and a non-stationary (with trend) Bayesian models for ﬂood frequency estimation were used, and their results were compared using the corresponding ﬂood frequency curves (FFCs), along with their uncertainty bounds. The results of trend analysis revealed signiﬁcant positive trends for 27.6% of the gauges, and 10% showed signiﬁcant negative trends at the signiﬁcance level of 0.05. In addition to these, 6.9% of the gauges also represented signiﬁcant positive trends at the signiﬁcance level of 0.1, while the remaining stations displayed insigniﬁcant trends. The non-stationary Bayesian model was found to be reliable for study sites possessing a statistically signiﬁcant trend at the signiﬁcance level of 0.05, while the stationary Bayesian model overestimated or underestimated the ﬂood hazard for these sites. Therefore, it is vital to consider the presence of non-stationarity for sustainable ﬂood management under a changing environment in the KRB, which has a rich history of ﬂooding. Furthermore, this study also states a regional shape parameter value of 0.26 for the KRB, which can be further used as an informed prior on shape parameter if the study site under consideration possesses the ﬂood type “ﬂash”. The synchronized appearance of a signiﬁcant increase and decrease of trends within very close gauge stations is worth paying attention to. The present study, which considers non-stationarity in the ﬂood regime, will provide a reference for hydrologists, water resource managers, planners, and decision makers.


Introduction
The comprehensive understanding of flood regimes is an important challenge in hydrology. Hydrologists and engineers customarily use flood frequency analysis (FFA) as a tool to understand area during 1999-2016 [28]. The Peshawar valley, with a rich history of flooding, provides the junctions for the Kabul River and its various right and left tributaries.
The above studies clearly show the presence of trends in rainfall regime as well as land use change in different sub-basins of the KRB. These climate and human interventions may induce non-stationarity in the flood regime. However, no studies have been reported to examine the presence or absence of stationarity in the flood regime of the basin. Therefore, it is imperative to study floods with a non-stationary point of view for the KRB.
Recently, Milly et al. [29] stated that the hypothesis of stationarity must be relinquished and that "stationarity is dead" and "should not be revived". The methods used for estimation of hydrologic indicators should be based on an innovative approach that would be reliable and useful for water management under a changing environment.
In the literature, various approaches have been reported using probabilistic modeling of flood frequency in a non-stationary context. Khaliq et al. [2] presented a comprehensive review, including the incorporation of trends in the parameters of the distributions, the incorporation of trends in statistical moments, the quantile regression method, and the local likelihood method. The studies of FFA under non-stationary conditions have mostly assumed trends in time [30][31][32][33][34][35][36][37]. The present study outlined a Bayesian framework for "at site flood frequency modeling" in stationary and non-stationary conditions. The fundamental concept is based on the generalized extreme value (GEV) distribution, combined with Bayesian inference for uncertainty assessment. For this study, a model with trend (non-stationary) and without trend (stationary) was used.
Previous studies in the KRB were limited to inundation mapping of flood-prone areas with a very little flow gauge station data, using a traditional frequentist approach [38][39][40][41][42].
The main objectives of the study were: (1) to analyze temporal and spatial trends in the annual maximum flood regime for the KRB, Pakistan, because no study has yet been reported in the literature to study the trends in annual extreme data of flood in detail, and (2) to address the non-stationary modeling of the flood regime in the KRB and its implications for flood management in a changing environment. We explored the differences between stationary and non-stationary flood quantile estimates for a given return period using flood frequency curves (FFCs), along with their uncertainty bounds for risk assessment, to analyze the importance of non-stationary models for improving flood management in the study area.

Study Area
The Kabul River basin (KRB), in Pakistan, stretches from 71 • 1 55"-72 • 56 0" E to 33 • 20 9"-36 • 50 0" N, as shown in Figure 1, which covers an area of 33,709 km 2 . The Kabul River starts at the base of Unai pass from the Hindu Kush Mountains in Afghanistan and flows eastward, covering a distance of 700 km to drain into the Indus River, Pakistan [43]. The entire basin covers an area of 87,499 km 2 . The elevation in the basin varies substantially from 249 m.a.s.l to 7603 m.a.s.l. High elevation mountains are mainly located in the north. The average temperature and average precipitation vary significantly across the River basin. The average temperature is about 13 • C. Most of the precipitation occurs in the northern mountain and highlands, reported up to 1600 mm. [44].
This study explores the part of the KRB that contributes to flooding. The flood problem arises mainly as the Kabul River enters Pakistan. The Logar River basin, Alingar River basin, and Panjshir River basin lie in Afghanistan. Three dams-Naghlu, Surobi, and Darunta-are located in Afghanistan on the Kabul River and Warsak dam is also located on the Kabul River in Pakistan. The study area is further divided into eight sub-basins: Kabul River basin, Chitral River basin, Main Swat River basin, Panjkora River basin, Lower Swat River basin, Kalpani River basin, Jindi River basin, and Bara River basin. The SRTM-DEM (Shuttle Radar Topography Mission-digital elevation model) of 30 m resolution and the geographical location of the sub-basins are also illustrated in Figure 1.

Flood Data
Twenty-nine flow gauge stations were selected to study the flood regime of the KRB. The annual maximum daily peak flow data for the seven flow gauge stations at main rivers sites were obtained from Surface Water Hydrology Project of Water and Power Development Authority (SWHP-WAPDA). The streamflow data of the remaining study sites were obtained from the Hydrology Section of the Irrigation Department of Khyber Pakhtunkhwa Province, Pakistan. The study sites that had at least 30 years of records were selected. The main characteristics of the sub-basins and the respective flow gauge stations in each sub-basin are presented in Table 1 and Figure 1 describes the geographical locations of the flow gauge stations in each sub-basin.

Flood Generating Mechanism in KRB
The hydrology of floods is linked to weather and climate as well as to physiographical features [45]. The basin has large altitudinal variations from 249 m.a.s.l. to 7603 m.a.s.l. Glacier-melt contribution from the upper part of the basin combined with rainfall in the lower part is the most likely cause of flooding in the region [38]. In the KRB, floods are mostly generated by monsoon rainfall but snow or glacial melt floods have also been observed in some parts of the basin. Snowmelt floods are not common. According to the data used in this study, all of the flood peaks were observed during the monsoon season, from July to August, in almost all the tributaries of the KRB. The historical floods occurred in July 2010, August 1995, and July 1992; all were observed during the monsoon. Anjum et al. [46] provided the details regarding rainfall magnitude, intensity, and spatial extent for the 2010 event. The South Asian monsoon originating from the Bay of Bengal is the dominant weather system for flood generation in the KRB.
However, the flood of 2005 in the Kabul and Indus Rivers was due to snowmelt as well as rainfall in the pre-monsoon period [47]. The flooding behavior of the different tributaries differs according to their catchment characteristics. The riverine floods in the Kabul River usually start below the Warsak dam, and this phenomenon propagates until its confluence with the Indus River at Khairabad

Flood Data
Twenty-nine flow gauge stations were selected to study the flood regime of the KRB. The annual maximum daily peak flow data for the seven flow gauge stations at main rivers sites were obtained from Surface Water Hydrology Project of Water and Power Development Authority (SWHP-WAPDA). The streamflow data of the remaining study sites were obtained from the Hydrology Section of the Irrigation Department of Khyber Pakhtunkhwa Province, Pakistan. The study sites that had at least 30 years of records were selected. The main characteristics of the sub-basins and the respective flow gauge stations in each sub-basin are presented in Table 1 and Figure 1 describes the geographical locations of the flow gauge stations in each sub-basin.

Flood Generating Mechanism in KRB
The hydrology of floods is linked to weather and climate as well as to physiographical features [45]. The basin has large altitudinal variations from 249 m.a.s.l. to 7603 m.a.s.l. Glacier-melt contribution from the upper part of the basin combined with rainfall in the lower part is the most likely cause of flooding in the region [38]. In the KRB, floods are mostly generated by monsoon rainfall but snow or glacial melt floods have also been observed in some parts of the basin. Snowmelt floods are not common. According to the data used in this study, all of the flood peaks were observed during the monsoon season, from July to August, in almost all the tributaries of the KRB. The historical floods occurred in July 2010, August 1995, and July 1992; all were observed during the monsoon. Anjum et al. [46] provided the details regarding rainfall magnitude, intensity, and spatial extent for the 2010 event. The South Asian monsoon originating from the Bay of Bengal is the dominant weather system for flood generation in the KRB.
However, the flood of 2005 in the Kabul and Indus Rivers was due to snowmelt as well as rainfall in the pre-monsoon period [47]. The flooding behavior of the different tributaries differs according to their catchment characteristics. The riverine floods in the Kabul River usually start below the Warsak dam, and this phenomenon propagates until its confluence with the Indus River at Khairabad near Attock. Riverine floods also occur in the Swat River in the Lower Swat catchment. In the rest of the KRB, flash flooding is a common disaster, along with landslides and torrential rains [45].

Trend Analysis
The non-parametric rank-based Mann-Kendall (MK) [48,49] test was used to detect trends in annual maximum flood series. The trend analysis was performed to show a clear understanding of the whole study area, while an objective criterion (means if statistically significant trend exists) was adopted for the non-stationary modeling of flood regime. The Mann-Kendall test was applied at different significance levels. The autocorrelation function (ACF) was also computed, before applying the Mann-Kendall (MK) test to check the presence of serial correlations in the annual maximum flood series.

Selection of Extreme Value Distribution
The Bayesian method using the GEV distribution is getting attention for analyzing hydrological extremes. The current study also utilized the GEV distribution, which is the integration of Gumbel, Fréchet, and Weibull distributions, and developed on the limit theorems for block maxima or annual maxima [50]. Mathematically, the cumulative distribution of the GEV can be written as [51]: where ψ(x) is expressed as 1 + ξ x−µ σ > 0; somewhere else, ψ(x) is either 0 or 1 [52]. The location parameter (µ), describes the center of the GEV distribution, the scale parameter (σ) describes the deviation around (µ), and the shape parameter (ξ) describes the tail behavior of the distribution. When ξ → 0, ξ < 0, and ξ > 0, GEV approaches the Gumbel, Weibull, and Fréchet distributions, respectively.

Goodness of Fit Statistics to GEV Distribution
The goodness of fit analysis of annual maximum peak flow data to the GEV distribution was performed in order to investigate whether the historical data belonged to the said GEV distribution. The Anderson-Darling (AD) [53] and Kolmogorov-Smirnov (K-S) [54] tests were performed for this purpose, using an EasyFit software (version 5.6, MathWave Technologies) [55]. EasyFit estimated the parameters of the GEV distribution based on maximum likelihood (ML) estimation, using equal probability sampling. The parameters estimated using the EasyFit software were used to assess the goodness of fit by AD and K-S statistics.
The K-S statistic (D) is based on the largest vertical difference between the theoretical and the empirical cumulative distribution function as shown below: The Anderson-Darling procedure compares the fit of an observed cumulative distribution function to an expected cumulative distribution function. This test gives more weight to tails than the K-S test.
H 0 : The data follow the specified distribution; H a : The data do not follow the specified distribution. The hypothesis regarding the distributional form was rejected at the significance level of 0.05 (alpha) if the test statistic, A 2 or D, was greater than the critical value of 2.5018 and 0.18482, respectively. Moreover, the outlier's detection in the annual extreme data of flood series was also performed using the Chauvenet's Criterion [56].

Model Design
The extreme value theory of stationary random process is based on that the statistical properties of extremes-here, the distribution parameters θ = (µ, σ, ξ) are free from time dependency [57], while in a non-stationary random process, the parameters of the said distribution function rely on time-dependency, and the properties of the distribution also vary with time [58]. For this study, two cases were considered.
(1) Stationary Case: all the model parameters were considered constant.
(2) Non-stationary Case: the location parameter (µ) was considered a function of time, as shown in Equation (4), while scale and shape parameters were kept constant: where t is time, θ = (µ 1 , µ 0 ) are the regression parameters [50,[57][58][59][60]. The location parameter was calculated for each study site in the stationary case and non-stationary case.

Bayes Theorem for GEV Distribution
Let θ be the parameter of given distribution and let Y = {y 1 , y 2 , . . . , y n } be the set of n observations. According to the Bayes theorem, the probability of θ given Y (posterior) is proportional to the product of the probability of θ (prior) and the probability of Y given θ (likelihood function). Assuming the independence between the observations, Y: Here, the likelihood function is the GEV distribution and θ is the vector containing the parameters of GEV distribution to be estimated. In the stationary case, θ = (µ, σ, ξ). By assuming independent GEV parameters: In the case of non-stationary analysis, θ contains an additional parameter, which is time-dependent here, i.e., µ(t), hence, the Bayes theorem for estimation of GEV parameters under the non-stationary case can be expressed as [57,60]: The resulting posterior distributions P(θ|Y, t) provide information on the distribution parameters (µ 1 , µ 0 , σ, ξ).

Prior Distribution
A Bayesian model utilizes a prior belief to calculate the posterior belief. For the current study, we utilized NEVA (non-stationary extreme value analysis, Matlab Package) [60][61][62] for our analysis. In NEVA, the priors are non-informative normal distributions, for location and scale parameters, while the priors for the shape parameter are a normal distribution, with a standard deviation of 0.3, as suggested by [57,60,62]. Initially, the shape parameter was considered a non-informative prior: if the value of the shape parameter in the posterior distribution exceeded beyond the plausible limit (−5, 5), as suggested by Martins and Stedinger [63]. Then, we modified the priors for shape parameter, considering partial pooling of information across sites that had similar flood types, for improving the flood quantiles estimates for "at site modeling" using the regional information. The shape parameter was considered an informative prior and the range of priors for shape parameter was: where, the Ksi stands for the shape parameter value of the site of interest from where it was exchanged. However, the location and scale parameter across sites were not shared.

Parameters Estimation and Convergence Criterion
To estimate the parameters inferred by Bayes, the Differential Evolution Markov Chain (DE-MC) is integrated to generate a large number of realizations from the parameters' posterior distributions [64,65]. The DE-MC attributes to the genetic algorithm Differential Evolution (DE) for global optimization over real parameter space with the Markov Chain Monte Carlo (MCMC) approach [64,65]. Here, the target posterior distributions were sampled through five Markov Chains constructed in parallel. These chains were allowed to learn from each other by generating candidate draws based on two random parent Markov Chains, rather than running independently. Therefore, it had the advantages of simplicity, speed of calculation, and convergence over the conventional MCMC. The initial numbers of burned samples were 6000 and numbers of evaluations were 10,000 for each study site. The R-hat criterion, suggested by Gelman and Shirley [66], was used to assess convergence, where R-hat should remain below 1.1.
Uncertainty estimates for FFCs are crucial for risk assessment and decision making. By combining DE-MC with Bayesian inference, the posterior probability intervals or credible intervals and uncertainty bounds of estimated return levels based on the sampled parameters could be obtained simultaneously for FFCs. For example, for a time series of annual maximum peak flow, the time-variant parameter (µ(t)) was derived by computing the 95th percentile of DE-MC sampled µ(t), (i.e., the 95th percentile of µ(t = 1), . . . , µ(t = 100)). These model parameters were then used to develop the stationary and non-stationary FFCs.
FFCs could also be drawn at 50% Bayesian credible intervals or at any other desired intervals.

Model Evaluation
In order to evaluate the suitability of the stationary versus non-stationary models, a Bayes factor K was calculated based on the posterior distributions of sampled parameters of both models. The stationary model was considered a null model M1, while the non-stationary model M2 was considered an alternative.
A value of Bayes factor > 1 denotes the stationary model is favored, while a value < 1 argues in the favor of the non-stationary model. Similarly, a value approaching +infinity favors the stationary model, and −infinity favors non-stationary models. Equation (10) represents the computation of Bayes factor, as follows: The term DA denotes input data, and θ stands for model parameters. The term Pr (DA|M) can be expressed using Monte Carlo integration estimation as follows: For more details see [67].

Temporal and Spatial Trends in Flood Regime
The trend magnitude for each station is presented in Table 2. Trend analysis results demonstrated the significant trend by 37.93% of the flow gauge stations in the entire basin; among them, 27.6% showed a significant increasing trend at the 0.05 significance level and 10.34% of the stations showed a significant decreasing trend at the significance level of 0.05. Moreover, 6.9% of the flow gauge stations also revealed a significant increasing trend at the significance level of 0.1. Non-significant trends were also exhibited by 31% of the flow gauge stations. The Chitral River at Chitral, the Kalpani River at Risalpur, the Kalpani River at Mardan, the Swat River at Chakdara, the Swat River at Ningolai, the Adezai River, the Naranji Nullah, the Bagiari Nullah, the Lund Khawar West, and the Bara River at Kohat Bridge displayed a significant increasing trend (Site #5, 6, 9, 11, 16, 18, 19, 21, 22, and 24), while the Swat River at Khawazakhela, the Naguman River, and the Badri Nullah showed a significant decreasing trend. The Khiyali River, the Panjkora River, and the Jundi Nullah at Tangi represented a moderate increasing trend, while the Kabul River at Warsak, the Swat River at Kalam, the Swat River at Munda Head Works, the Budni Nullah, and the Khuderzai Nullah displayed a moderate decreasing trend. Moreover, the main flow gauge station-the Kabul River at Nowshera-did not show any significant trend. *** Trend is significant at α = 0.001, ** Trend is significant at α = 0.01, * Trend is significant at α = 0.05, + Trend significant at α = 0.1. Figure 2 represents the basin-wide spatial distribution of trends in the flood regime, which showed that the flood regime of the Chitral River, the Kalpani River, and the Main Swat River basins exhibited significant increasing trends. However, the Swat River at Khawazakhela and the Badri Nullah represented a significant decreasing trend.
The Lower Swat River, the Kabul sub-basin, and the Jindi River basins showed non-significant trends. Especially for the southwestern part of the KRB, the Bara River at Kohat Bridge showed non-stationarity in the flood regime due to a significant increasing trend, while all other flow gauge stations in the Bara River basin showed insignificant trends.
The change in flood regime was found to be more evident for the northern and northeastern part of the KRB as compared to the central and southwestern parts of the KRB. Consequently, the overall basin showed large spatial variations. However, the basin did not represent a regular spatial pattern. These spatial variations may be due to climatology, topography, and complex orography of the KRB in the HKH region. However, the temporal changes in the flood regime might be attributable to regional environmental change. The results of trend analysis were consistent with previous studies [68][69][70], but these studies used only two to three flow gauges.

Evaluation of Goodness of Fit for Annual Extreme Data of Flood
An objective criterion as suggested by Rosner et al. [71] was adopted to evaluate the goodness of fit of the annual maximum flood series, with GEV and non-stationary temporal trend modeling of the flood regime, i.e., only the sites showing significant trends in their flood regime were selected. The study sites under consideration showed proper fitting using the AD test and the value of test statistics for all the sites using the AD test was less than the critical value of 2.5018. While using the K-S test, all the sites showed proper fitting except site 21, but it was included in the analysis because it satisfied the AD test. Table 3 provides the results of the test statistics and estimated GEV parameters using ML. The p-value belongs to K-S test only.
Outliers were also detected in the data of the annual maximum flood series for the KRB. Table 4 demonstrates the outliers present in the data of selected study sites. Sites 5,6,16,18,19,21,

Evaluation of Goodness of Fit for Annual Extreme Data of Flood
An objective criterion as suggested by Rosner et al. [71] was adopted to evaluate the goodness of fit of the annual maximum flood series, with GEV and non-stationary temporal trend modeling of the flood regime, i.e., only the sites showing significant trends in their flood regime were selected. The study sites under consideration showed proper fitting using the AD test and the value of test statistics for all the sites using the AD test was less than the critical value of 2.5018. While using the K-S test, all the sites showed proper fitting except site 21, but it was included in the analysis because it satisfied the AD test. Table 3 provides the results of the test statistics and estimated GEV parameters using ML. The p-value belongs to K-S test only.
Outliers were also detected in the data of the annual maximum flood series for the KRB. Table 4 demonstrates the outliers present in the data of selected study sites. was also recorded as an outlier. Despite the existence of outliers, the data series belongs to the GEV distribution as per AD and K-S test statistics. Table 3. The goodness of fit statistics of annual maximum daily peak flow to generalized extreme value (GEV) distribution.

Regionalization of Shape Parameter for Flash Floods Across the KRB
The shape parameter of the GEV distribution is important for the estimation of flood quantiles. Initially, non-informative priors for shape parameter were used for the Bayesian analysis of annual maximum flood regime for the selected study sites that had non-stationarity, due to the existence of temporal trends. About 30% of the study sites yielded a value of shape parameter in posterior distribution that exceeded beyond the plausible limit (−5, 5) as suggested by Martins and Stedinger [63], causing the degeneration of the GEV distribution. In order to avoid this, homogeneous sites were  [72][73][74][75][76] state that the use of regional information will improve the flood frequency estimation and reduce the uncertainty for sites having short records. Table 5 illustrates the correlation matrix for the selected study sites, which demonstrates that hierarchical clustering is possible based on the correlation between the annual maximum flood series of the selected study sites. A positive correlation is present between the sites. Although site 9 possesses the least positive correlation with site 6, site 6 has more positive correlations with other sites, hence why all the sites are considered homogenous. All the sites also possess "flash" flood type. Moreover, all the sites could also be considered homogenous because of the existence of trends in their flood regime. Sun et al. [77] also highlighted the clustering of temporal trends and exchange of shape parameter for the Bayesian analysis of annual maximum floods across Germany.
Furthermore, the utilization of ML for the estimation of the shape parameter for GEV distribution was found reliable for large records-at least 50 year [78]. After considering all the study sites as homogenous based on the correlations between sites, similar flood type, and existence of trends, Naranji Nullah (site 16), with a sufficiently long record, was considered the benchmark site. The shape parameter estimated by ML was approximately 0.26 for site #16. This value of shape parameter (0.26) was further recognized for all the study sites as an informative prior in the Bayesian model. This is like partial pooling of information across homogenous sites, which ultimately improved the flood quantiles estimates using the regional information as compared to non-informative priors on shape parameter. Lima et al. [79] used the basin's average shape parameter value in local and regional hierarchical Bayesian models to solve the issue of sites where the shape parameter value exceeds beyond (−5-5). Lima et al. [79] used the prior for shape parameter as non-informative, but this study considers the priors on shape parameter to be informative priors.

Comparison between Stationary and Non-Stationary Bayesian Models
Non-stationary FFCs were constructed for the flow gauges with significant increasing trends in their flood series and compared with their stationary FFCs, considering the entire data series without the elimination of outliers. Table 6 demonstrates the results at the 95% Bayesian credible interval for all the selected study sites in the KRB. The stationary model showed overestimation as compared to the non-stationary model for a 100-year flood (The flood having a probability of exceedance of 0.01) by 1494 m 3 s −1 (+34.9%) for the Adezai River at Adezai Bridge (site 5). The value of the Bayes factor was 0.0058, which was less than 1, so the non-stationary model was favored. The maximum peak flood observed at the Adezai River was 2285 m 3 s −1 , during the historic flood of 2010 in the KRB. The non-stationary model reasonably described the historical peak flood (Figure 3). Similar behavior was observed for the other study sites, 11, 18, 21, 22, and 24 (Table 6 and Figures 4-8).               On the other hand, the stationary model also underestimated the 100-year flood as compared to the non-stationary model. For example, the stationary model underestimated the 100-year flood by 23 m 3 s −1 (−1.19%) as compared to the non-stationary model ( Figure 9, the Chitral River at Chitral).
The value of the Bayes factor was 0.068, which was less than 1, so the non-stationary model was favored. This behavior was obvious for study sites 6 and 19 (Table 6, Figures 9 and 10).     Furthermore, for study sites possessing trends at the significance level of 0.1 (although the trends were modeled at 0.1% for these sites), the non-stationary model overestimated the 100-year flood as compared to the stationary model, and the corresponding value of the Bayes factor was greater than 1. This ultimately favors the stationary model for sites 9 and 16 (Table 6, Figures 11 and 12).  The non-stationary model was found to be reliable for the sites to study the annual maximum flood regime, which exhibits significant trends at α = 0.05. The uncertainty bounds of most of the  The non-stationary model was found to be reliable for the sites to study the annual maximum flood regime, which exhibits significant trends at α = 0.05. The uncertainty bounds of most of the The non-stationary model was found to be reliable for the sites to study the annual maximum flood regime, which exhibits significant trends at α = 0.05. The uncertainty bounds of most of the study sites were high because of the higher value of coefficient of variation, except site 6, where the value of coefficient of variation was 0.2 (see Table 1 and Figure 9). The uncertainty bound could also be higher due to the existence of outliers in the data series. The removal of outliers in the data series can reduce uncertainty.

Performance of Bayesian Models to Predict the Extreme Floods
The Bayesian models were again re-run for all the study sites before the commencement year of the extreme flood event that we labeled an outlier for all the study sites. Similarly, the objective criterion was adopted again (if a trend exists then modeling was performed using the non-stationary Bayesian model, otherwise only the stationary Bayesian model was used). Table 7 describes the results of the corresponding stationary and non-stationary Bayesian models.
The Bayesian models performed well, predicting the extreme floods satisfactorily for all the study sites except sites 9 and 22 (see Table 7 and Figures 13-22). The stationary Bayesian model was found reliable for sites 5, 6, 11, 16, 21, and 24. However, the non-stationary model was favored as compared to the stationary model as per the Bayes factor criterion (see Table 7) for sites 6, 18, 19, and 21. Despite the existence of a significant trend at site 11, the stationary model was favored as per the Bayes factor criterion, as well as predicting the extreme flood of 2016.
The outlier's removal improved the fitting of the FFCs as compared to considering the entire data series and also reduced the uncertainty (see . Statistical modeling of extremes in hydrology is exciting and challenging, and opens the door for further studies. For example, for study site 9 (the Swat River at Chakdara), it is better to consider the entire data series, and for better fitting, the incorporation of monsoon rainfall as a covariate might be fruitful. The modeling could also be performed by considering other distributions in the Bayesian framework or using the traditional frequentist framework.
Finally, from a regional perspective, the region is heterogeneous due to large altitudinal variations. Due to the regional heterogeneity associated with elevation, it seems to be quite difficult to develop a regional Bayesian model for the whole KRB, but efforts should be made to develop a regional Bayesian model at sub-basin or catchment scales in future studies by further pooling of information for other parameters, like location and scale, across sites. Moreover, the studies are also required to deeply understand the impact of climate or dominating weather patterns, such as the South Asian monsoon, low climate variability El-Niño Southern Oscillation and Indian Ocean Dipole (ENSO, IOD, etc.) and human factors, such as land use cover change (LUCC), population increase, and reservoir construction, on the flood regime of the KRB [1,13,[80][81][82][83].                    Statistical modeling of extremes in hydrology is exciting and challenging, and opens the door for further studies. For example, for study site 9 (the Swat River at Chakdara), it is better to consider the entire data series, and for better fitting, the incorporation of monsoon rainfall as a covariate might   Statistical modeling of extremes in hydrology is exciting and challenging, and opens the door for further studies. For example, for study site 9 (the Swat River at Chakdara), it is better to consider the entire data series, and for better fitting, the incorporation of monsoon rainfall as a covariate might

Conclusions
Analyzing the flood regime and its non-stationary modeling in a Bayesian framework for the KRB was the main objective of the present study. To achieve this, a Mann-Kendall trend analysis was performed to explore the flood regime of the KRB in detail, and finally, the stationary and non-stationary Bayesian models with informed priors on shape parameter for GEV distribution were developed, and their results were compared by using the corresponding FFCs, along with their uncertainty bounds. We utilized the annual extreme data of flood series for the study area, with a maximum record of 1962-2016 (55 years) and a minimum of 1987-2016 (30 years). The key findings of the study are described below:

1.
Trend analysis showed a mixture of increasing and decreasing trends at different gauges in the KRB at α = 0.05. The Chitral River, Kalpani River, Main Swat River, and Bara River basins showed significant increasing trends, and the Panjkora River basin displayed a moderate increasing trend in its annual maximum flood regime. However, the Lower Swat and Kabul sub-basins showed decreasing trends, except for the Adezai River in the Kabul sub-basin, which showed a significant increasing trend. 2.
The overall basin was under critical change and signals of clear non-stationarity in the flood regime were evident at various spatial scales throughout the basin.

3.
The presence of a significant trend and significant difference in flood estimates for 100-year flood between stationary and non-stationary FFCs were found that represent the clear violation from the so-called stationary assumption.

4.
The non-stationary Bayesian model was found to be reliable for the study sites that had a significant trend at α = 0.05, while the stationary model overestimated or underestimated the flood risk for these sites. On the other hand, the stationary Bayesian model performed better for the study sites for trends at α = 0.1, while the non-stationary Bayesian model overestimated or underestimated the flood risk for such sites. 5.
The use of informed priors on the shape parameter based on regional information improved the estimation of flood quantiles and reduced the uncertainty. 6.
Proper consideration should be given to identify the outliers while using Bayesian models. 7.
The presence of non-stationarity in the flood regime of the KRB has substantial implications for flood management and water resources development. A design with stationary assumption will cause two major concerns: under estimation or overestimation of design for structural and non-structural measures in the KRB. An event-based design may also overestimate or underestimate the risk in hydraulic design that was intended. Some previous studies in other parts of world also provided similar results [1,13,31,[84][85][86].
The study will be helpful for sustainable flood management and provide a reference for studying floods in a changing environment for hydrologists, water resources managers, decision makers, and concerned organizations.