Introduction

Floods are the most destructive natural catastrophes, leading to losses in human lives, infrastructure, crop yields, and industrial production. It can also seriously affect the environment and the social lives of people in terms of electrical damages and risks, contaminated water, sanitary hazards, problems in drainage, and destruction of roads and other communication networks. A flood is a situation when water overflows from a watercourse (river, lake, or stream channel) and comes to the land, which is usually dry [1-4]. The main reasons behind floods are heavy rains, melting glaciers, breakage in reservoirs used for water storage, and a lack of river channels to effectively convey the excess water [5-7]. Floods are natural phenomena occurring at irregular intervals of time. The flood frequency analysis is used with the objective of predicting the recurrence of flood levels for different return periods. Such prediction of floods is important for certain planning at the government and country level, like the construction of water reservoirs, spillways, bridges, canals, and culverts. It is also important to plan disaster management activities (including protection of human lives, assets, and food security) for any possible hazard in the future. Flood frequency can be carried out with different approaches, like the California method, the Hazen method, and Gumbel’s method [1]. The flood analysis can be carried out at the country, regional, or site level. At-site flood frequency analysis is the most appropriate, effective, and simple way to predict flood recurrence at a particular gauging station [5, 6, 8-11].

The choice of the probability model blended with an appropriate estimation strategy plays a vital role in the modeling of flood data. Each distribution has its own advantages and disadvantages in describing flood data. There can be different probability distributions with different parameter estimation methods suitable for modeling maximum water flow data from different gauging sites. In the literature concerning at-site flood frequency analysis, different probability distributions like Generalized Extreme Value (GEV), GEV type-1, GEV type-2, normal, generalized normal, two-parameter lognormal, three-parameter lognormal, two-parameter gamma, three-parameter gamma, Pearson type-III, Gumbel, reverse Gumbel, exponential, generalized logistic, four-parameter Wakeby, five-parameter Wakeby, and generalized Pareto have been applied for flood frequency analysis of different rivers in different countries [1, 5, 6, 12-18]. A brief review can be seen in some studies [19-21].

Following the findings of many studies for different rivers in the world, we have used Generalized Logistic (GLO), Generalized Extreme Value (GEV), Gumbel (GUM), Generalized Pareto (GPD), and Reverse Gumbel (REV.GUM) distributions for at-site modeling of annual maximum flood data at 5 sites of the Chenab River. Some of these probability models were also used by Ahmad et al. [22] for at-site analysis of many gauging sites, including those of the Chenab River, but with relatively smaller sequences of peak flow data.

Several estimation methods are available for parameter estimation of probability models used for flood frequency analysis. In the present study, we have applied L-Moments (LM), Maximum Likelihood (ML), Least Squares (LS), Weighted Least Squares (WLS), and Relative Least Squares (RLS) methods for estimating the parameters of considered probability distributions.

The performance of each probability distribution estimated with a particular estimation method is assessed based on different goodness-of-fit tests and accuracy measures.

The present study is focused on at-site flood frequency analysis of the Chenab River, with the prime objective of identifying the most suitable probability distribution with a parameter estimation method for modeling the annual maximum discharge data of different gauging sites of the Chenab River. After identifying the most appropriate probability distribution and estimation method for each gauging station, the maximum annual peak flow for different return periods with a specific non-exceedance probability has been estimated.

The rest of the article is structured as follows: Section 2 provides a brief description of the Chenab River and the data used in the current work. Section 3 describes the candidate probability models and parameter estimation methods. Section 4 presents and discusses the results obtained by applying different probability distributions with estimation methods at different gauging sites. Finally, conclusions drawn from the study are presented in Section 5.

Data and Gauging Sites

The Chenab River is one of the major rivers in Punjab Province, Pakistan. The river flows from northeast to southeast directions and enters Pakistan’s Punjab province from Indian Punjab. It measures 960 Km with an approximate enactment area of 29000 km2. Its average annual flow is 12.38 MAF [23]. A map of the Chenab River showing the considered gauging sites is reproduced from Magsi and Atif [24] in Fig. 1. The data have been collected from the Hydrological Directorate of Discharge and Flood Zone, Lahore. The range of the data is different for different gauging sites, varying from 48 to 90 years. The longitude, latitude, elevation, and length of the data series for each gauging site are given in Table 1. The maximum annual flow at each gauging site is graphically represented in the form of time plots in Fig. 2.

Materials and Methods

The choice of a probability distribution and estimation method is vital in any at-site flood frequency analysis. In the following sub-section, a brief introduction of candidate probability distributions and estimation methods used in the current study is provided.

Candidate Probability Distributions

The choice of a suitable probability distribution is vital in any at-site flood frequency analysis. Many probability distributions have been used for flood frequency analysis in different countries. For the analysis of five gauging sites of the Chenab River, we have used Generalized Logistic (GLO), Generalized Extreme Value (GEV), Gumbel (GUM), Generalized Pareto (GPD), and reverse Gumbel (REV.GUM) probability distributions. These distributions are applied and found plausible for at-site flood frequency analysis of different rivers and also for other environmental variables [14, 25-30].

L-Moments (LM) Method

Introduced by Hosking [31, 32], the L-moments are linear combinations of probability-weighted moments (PWM). Being the direct measure of the scale and shape of the data, L-moments are more suitable and convenient than probability-weighted moments. However, it is common practice to compute L-moments using PWM [5, 33, 34]. Theoretically, the rth PWM βr, for a probability distribution with a quantile function φ(F), is defined as:

βr=∫01ϕ(F) FrdFr=0, 1, 2,...
Fig. 2.

Time plots of maximum annual discharge at different gauging sites.

https://www.pjoes.com/f/fulltexts/185970/PolJEnvironStud-2025-34-1-29-42_Fig_2_min.jpg
Fig. 1.

A map showing gauging sites of the Chenab Rive (Source: Magsi and Atif, 2012).

https://www.pjoes.com/f/fulltexts/185970/PolJEnvironStud-2025-34-1-29-42_Fig_1_min.jpg
Table 1.

Brief of gauging sites.

Gauging siteTime periodLongitudeLatitudeElevation (meter)
Marala1946-201774.464332.6723234
Khanki1944-201773.968732.3988219
Qadirabad1970-201773.686332.3193205
Trimmu1928-201772.146231.1447150
Punjnad1948-201771.019829.34641402

For sample data, the first four PWMs are computed as:

β0=n−1∑i=1nxi:n,β1=n−1∑i=2n(i−1)(n−1)xi:n,
β2=n−1∑i=3n(i−1)(i−2)(n−1)(n−2)xi:n,
β3=n−1∑i=4n(i−1)(i−2)(i−3)(n−1)(n−2)(n−3)xi:n.

Consequently, the first four L-moments are computed by the following relationships,

λ1=β0, λ2=2β1−β0, λ3=6β2−6β1+β0,
λ4=20β3−30β2+12β1−β0.

Hosking [32] defined the coefficients of variation (τ2), skewness (τ3) and kurtosis (τ4) in terms of ratios of the L-moments as:

τ2=λ2λ1,τ3=λ3λ2,τ4=λ4λ2.

Just like the method of moments, parameter estimation through the L-moments approach is done by solving a system of simultaneous equations obtained by comparing theoretical and corresponding sample L-moments. The LM method has gained much attention recently due to its computational simplicity [5, 6, 12, 35].

Maximum Likelihood (ML) Method

Maximum Likelihood (ML) is a widely applied estimation method. It gives those values as parameter estimates for which the likelihood, or log-likelihood, attains its maximum. Let us have a random sample of n observations x1, x2,... xn from probability density function f(xi|θ) where θ is the vector of unknown parameters. Let the likelihood and log-likelihood of n independent observations are L(θ)=∏i=1nf(xi|θ) and l(θ)=∑i=1nlogf(xi|θ), respectively. The ML estimate of θ is the value for which l(θ) attains its maximum. The maximization is achieved by equating the derivative of l(θ) with respect to θ to zero. If this process does not yield a closed-form solution, then maximization can be achieved numerically through an optimization technique. For the current work, we have applied algorithms available in R-language [36] for numerical optimization.

Least Squares Method

The Least Squares (LS) method is commonly used for parameter estimation in probability models [37- 39]. It is based on minimizing the following objective function:

LS(θ)=∑i=1n[F(xi)−in+1]2,

where F(xi) is the theoretical CDF of the corresponding probability distribution while i/(n+1) is the sample or observed CDF of the data. This can be done by solving the system of equations obtained by taking derivatives of the above function with respect to unknown parameters. The details can be seen in Ali [37]. For complex functions, the parameter estimates can be obtained by a suitable optimization algorithm.

Weighted Least Squares Method

As proposed by Bergman [40], Weighted Least Squares (WLS) is a variant of the LS method. Instead of identical weights for all observations, WLS uses different weights for different observations. Specifically, in the WLS method, parameter estimates are obtained by minimizing the following function:

WLS(θ)=∑i=1nWi[F(xi)−in+1]2.

Following different studies [41-44], we have used weight function as:

Wi=(n+1)2(n+2)i(n−i+1).

Relative Least Squares Method

Proposed by Pablo and Bruce [45], Relative Least Squares (RLS) is another variant of the LS method. Instead of squared errors (as in the LS method), the RLS method is based on minimizing the squared relative errors. Recently, RLS is reported to be a better performer, in the case of some two-parameter probability distributions [46-48]. Parameter estimates through the RLS method can be obtained by minimizing the following function:

RLS(θ)=∑i=1n[F(xi)−in+1in+1]2.

Performance Indicators

To suggest the most suitable probability model for a particular gauging site, the performance of candidate models with an estimation method is assessed. The performance comparison is carried out using some common measures of accuracy and goodness-of-fit (GOF) tests. The test statistics of the GOF tests used in the current study are defined as:

Kolmogorov-Smirnov (KS) test,

KS=Max{|Fn(xi)−F^(xi)|}·

Anderson-Darling (AD) test,

CVM=112n+∑in(2i−12n−F^(xi))2.

Cramér–von Mises (CVM),

AD=−n−∑i=1n(2i−1n{log(1−F^(xn−i+1))+log(F^(xi))}).

Similarly, accuracy measures like Root Mean Square Error (RMSE), Root Mean Square Percentage Error (RMSPE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE) are defined as:

RMSE=∑i=1n(xi−x^i)2n,RMSPE=∑i=1n(xi−x^ixi)2*100n,
MAE=100n∑i=1n|xi−x^i|,MAPE=100n∑i=1n|xi−x^ixi|.

In the above GOF tests and accuracy measures, Fn(xi) represents the observed or empirical CDF and F̂ (xi) is the expected CDF of the distribution with particular parameter estimates. Similarly, xi represents the observed value, while x̂i is the estimated value using particular parameter estimates.

Estimation of Quantiles

The estimation of maximum peak flow for a particular return period (T) xT is a major objective of any flood frequency analysis. The recurrence interval T is expressed in terms of years, in which the maximum peak flow xT is expected to return. A specific flood level xT (T being the return period) may be exceeded once in T years. Therefore, P(X>xT)=1T, then the cumulative probability of non-exceedance F(xT) is

F=F(xT)=P(X≤xT)=1−P(X>xT)=1−1/T

The quantile xT can be estimated by putting F in the quantile function of the corresponding probability distribution. If the quantile function of a probability distribution cannot be expressed in closed form, then the numerical relationships between xT and F can be used to evaluate xT for a given value of F. The quantile estimate for a return period of T years is just computed by inserting value F = (1–1/T) in the quantile function. The standard errors of these flood levels for different return years were estimated, and hence 95% confidence intervals were constructed using the parametric bootstrapping approach [5, 6, 49].

Results and Discussion

The descriptive statistics of maximum annual discharge data from different gauging stations are given in Table 2. From these summary measures, it can be seen that, on average, Khanki gauging station has the maximum annual peak flow, and this particular site also has a higher variation in the maximum yearly flow. Generally, we see that variability in annual maximum peak flow is not too different for different gauging stations, which is evident from values of the coefficient of variation ranging from 66.45% to 70.25%. Randomness, independence, stationarity, and skewness are the core assumptions for the data to be used in flood frequency analysis through probability models [6, 49, 50]. In the current article, we used the Wald-Wolfowitz (WW) test for testing randomness and the correlation coefficient at lag-1 for assessing independence, while the assumption of stationarity was tested using the Augmented Dickey-Fuller (ADF) test. The results of testing these assumptions are given in Table 3. From the p-values of the ADF and WW tests in Table 3, it is clear that assumptions of stationarity and randomness seem reasonable for all gauging sites considered in the study. Moreover, values of serial correlation at lag-1 do indicate that there is no serious/severe violation of the independence assumption. Similarly, values of the coefficient of skewness indicate that the data series for all gauging stations is positively skewed.

The validity of these assumptions justified the application of considered probability distributions for modeling of the annual maximum peak flow data for all gauging sites.

The parameter estimates through LM, ML, LS, WLS, and RLS methods for all of the candidate probability distributions are given in Tables 4-8 for gauging sites at Marala, Khanki, Trimmu, Qadirabad, and Punjnad, respectively. The parameter estimates are accompanied by standard errors (in parenthesis) computed using the parametric bootstrapping approach [5, 6, 49]. The comparative performance of different probability distributions with estimation methods is also given in Tables 4-8. The performance is assessed in terms of the total number of ranks based on three GOF tests and four accuracy measures. For each gauging site, the blend of probability model and estimation method is compared by ranking based on GOF tests and accuracy measures. The highest rank is given to the highest p-value (in the case of GOF tests) and the lowest value (in the case of accuracy measures). The total rank score computed by adding all ranks is given for different combinations of probability distribution and estimation methods. The probability distribution with the maximum total rank is considered to be the most suitable for the particular gauging site.

Table 2.

Descriptive statistics of annual maximum discharge (m3/s) of the Chenab River.

Gauging sitenMeanMedianMaximum (Annual)SDCV
Marala729418.877080.2931148.536315.720.6705
Khanki7410569.088004.8530765.127424.820.7025
Qadirabad4810231.018417.4526859.107170.630.7009
Trimmu908452.827326.7726731.105671.740.6709
Punjnad708208.127551.8622724.725454.290.6645
Table 3.

Results of testing assumptions of data from different gauging sites.

Gauging sitenADF(P-value)WW(P-value)Correlation at lag-1Skewness
Marala72<0.010.23530.0011.2856
Khanki74<0.010.81490.0471.1508
Qadirabad48<0.011.000.0790.9622
Trimmu90<0.010.39640.2871.2185
Punjnad70<0.010.05410.2950.6158

From these tables, it turns out that Generalized Pareto (GP) distribution with LS estimation is more suitable for Marala and Khanki gauging stations, as this blend of probability distribution and estimation method got the highest rank based on seven different performance indicators. For gauging sites at Qadirabad and Punjnad, GP distribution with the WLS method provides more accurate and precise results. Generalized extreme value distribution with ML estimation is best suited for the gauging site at Trimmu. It is noted that the GP distribution with LM and WLS estimations performs equally well for Punjnad.

Table 4.

Parameter estimates and performance comparison (Gauging site: Marala).

DistributionMethodμ̂α̂k̂Total Rank
GLOLM7767.663 (612.488)2811.449 (331.925)-0.316 (0.092)102
MLE7685.134 (5.892)3105.981 (5.67)-0.55 (0.035)83
LS7298.628 (8.7)3038.087 (8.322)-0.561 (0.138)113
WLS7324.065 (5.009)2902.162 (4.628)-0.503 (0.076)136
RLS7355.782 (4.984)3140.995 (4.89)-0.582 (0.034)96
GEVLM6230.542 (529.408)3770.856 (479.633)-0.216 (0.111)126
MLE6231.38 (9.024)3769.103 (8.699)-0.398 (0.12)107
LS6210.904 (10.72)3787.61 (9.787)-0.348 (0.231)120
WLS6222.21 (7.434)3777.271 (6.179)-0.33 (0.178)122
RLS6230.72 (9.075)3773.934 (7.98)-0.488 (0.188)72
GUMLM6642.637 (616.005)4809.696 (504.781)---73
MLE6645.577 (10.808)4806.611 (11.215)---70
LS6629.655 (12.357)4803.194 (10.47)---79
WLS6643.661 (14.906)4791.505 (12.161)---76
RLS6652.114 (11.902)4809.471 (11.106)---67
GPDLM2622.268 (296.351)7059.481 (1346.109)0.039 (0.137)142.5
MLE2637.917 (6.597)7054.279 (7.41)0.067 (0.086)136.5
LS2608.045 (11.247)7061.501 (10.225)0.014 (0.173)146.5
WLS2614.224 (6.396)7062.709 (5.011)0.013 (0.129)144.5
RLS2619.634 (10.667)7061.501 (9.184)-0.059 (1.293)125.5
REVLM12195.101 (641.604)4809.696 (511.94)---22.5
MLE12191.336 (8.003)4802.985 (8.62)---28.5
LS12190.788 (11.637)4808.515 (10.064)---30.5
WLS12190.117 (15.513)4823.409 (11.786)---27.5
RLS12194.263 (10.3)4798.256 (12.029)---28.5
Table 5.

Parameter estimates and performance comparison (Gauging site: Khanki).

DistributionMethodμ̂α̂k̂Total Rank
GLOLM8578.07 (667.167)3297.339 (386.536)-0.323 (0.093)96
MLE8327.079 (4.872)3719.802 (4.996)-0.641 (0.027)77
LS7829.055 (8.756)3603.425 (8.066)-0.65 (0.135)97
WLS7824.397 (4.275)3305.906 (4.286)-0.569 (0.069)116
RLS7804.565 (4.714)3520.225 (4.621)-0.649 (0.029)110
GEVLM6781.488 (596.143)4398.43 (528.019)-0.226 (0.11)115
MLE6779.137 (8.82)4387.121 (8.658)-0.554 (0.115)79
LS6770.721 (10.727)4406.267 (9.219)-0.37 (0.247)114
WLS6782.714 (7.262)4397.79 (6.158)-0.355 (0.175)115
RLS6760.904 (9.421)4379.058 (7.609)-0.634 (0.17)59
GUMLM7287.04 (708.425)5685.983 (611.485)---73
MLE7286.831 (11.241)5671.209 (11.493)---81
LS7281.227 (12.665)5693.248 (11.419)---77
WLS7287.752 (14.762)5669.175 (12.005)---77
RLS7289.501 (11.638)5685.56 (11.303)---67
GPDLM2597.849 (372.489)8150.795 (1491.281)0.022 (0.125)164
MLE2622.528 (6.373)8140.975 (7.507)0.065 (0.08)140
LS2606.5 (10.996)8156.522 (10.048)0.005 (0.17)170
WLS2597.743 (6.112)8152.081 (4.5)-0.001 (0.133)160
RLS2598.673 (11.097)8154.775 (8.818)-0.012 (0.502)153
REVLM13851.117 (653.604)5685.983 (577.742)---24
MLE13847.455 (8.518)5680.258 (9.077)---28
LS13840.106 (12.046)5692.337 (9.851)---27
WLS13812.404 (14.578)5687.364 (12.066)---35
RLS13855.784 (11.213)5682.703 (12.367)---21
Table 6.

Parameter estimates and performance comparison (Gauging site: Qadirabad).

DistributionMethodμ̂α̂k̂Total Rank
GLOLM8532.595 (943.523)3447.037 (482.036)-0.274 (0.101)89
MLE8140.054 (6.105)3823.902 (6.241)-0.592 (0.043)81
LS7916.106 (8.888)3906.402 (8.014)-0.544 (0.165)103
WLS7880.047 (5.242)3581.542 (5.366)-0.488 (0.102)119
RLS7824.25 (5.756)4133.924 (5.76)-0.667 (0.04)64
GEVLM6609.171 (760.99)4781.896 (676.025)-0.156 (0.122)112
MLE6617.385 (9.528)4785.468 (9.1)-0.372 (0.145)97
LS6602.378 (10.34)4783.172 (8.941)-0.283 (0.268)116
WLS6589.003 (7.828)4786.033 (6.999)-0.263 (0.22)126
RLS6603.9 (9.475)4783.234 (8.642)-0.508 (0.242)68
GUMLM6973.139 (837.479)5644.117 (760.89)---70
MLE6960.049 (11.416)5634.125 (11.208)---78
LS6954.55 (11.025)5636.569 (10.718)---79
WLS6984.612 (13.439)5635.43 (11.65)---70
RLS6972.529 (11.505)5632.61 (11.148)---78
GPDLM1858.87 (518.891)9544.299 (2122.91)0.14 (0.154)162
MLE1866.996 (6.332)9536.78 (7.265)0.206 (0.084)142
LS1841.839 (10.652)9522.373 (9.718)0.104 (0.179)160
WLS1864.111 (6.614)9538.813 (5.112)0.122 (0.141)171
RLS1860.427 (11.3)9535.66 (9.081)0.112 (1.109)165
REVLM13488.884 (876.851)5644.117 (709.547)---20
MLE13488.048 (8.474)5643.668 (9.037)---25
LS13489.536 (12.453)5659.823 (9.925)---26
WLS13468.959 (13.579)5629.302 (11.608)---32
RLS13508.622 (10.374)5615.216 (12.053)---22
Table 7.

Parameter estimates and performance comparison (Gauging site: Trimmu).

DistributionMethodμ̂α̂k̂Total Rank
GLOLM7137.086 (520.59)2647.512 (267.601)-0.276 (0.075)112
MLE6884.571 (6.052)2668.764 (6.829)-0.388 (0.041)129
LS6859.404 (9.089)2769.048 (8.109)-0.419 (0.123)118
WLS6903.337 (4.801)2732.891 (4.616)-0.375 (0.072)137
RLS6985.346 (5.814)2810.131 (5.347)-0.39 (0.046)125
GEVLM5661.177 (455.562)3666.866 (382.674)-0.158 (0.093)142
MLE5658.47 (8.095)3650.411 (8.223)-0.205 (0.09)149
LS5646.884 (11.211)3656.956 (8.9)-0.232 (0.2)129
WLS5656.272 (7.13)3666.872 (5.656)-0.217 (0.143)136
RLS5650.45 (9.273)3657.876 (8.406)-0.226 (0.164)135
GUMLM5945.849 (497.13)4343.213 (432.186)---67
MLE5945.236 (11.123)4335.3 (11.461)---73
LS5928.813 (12.522)4354.189 (10.629)---76
WLS5941.38 (17.229)4334.401 (13.378)---78
RLS5958.096 (11.527)4330.818 (10.811)---66
GPDLM2025.099 (273.118)7296.174 (1201.291)0.135 (0.114)106.5
MLE2027.821 (5.886)7307.77 (6.828)0.16 (0.064)97.5
LS2016.544 (11.502)7286.063 (10.02)0.128 (0.138)108.5
WLS2017.511 (5.55)7294.157 (4.441)0.116 (0.108)108.5
RLS2013.246 (10.849)7300.421 (9.2)0.281 (0.678)64.5
REVLM10959.79 (511.152)4343.213 (408.4)---22.5
MLE10959.79 (8.378)4343.213 (9.275)---22.5
LS10961.479 (11.472)4346.723 (10.155)---18.5
WLS10948.152 (17.279)4340.959 (13.053)---30.5
RLS10963.743 (10.636)4330.726 (11.514)---23.5
Table 8.

Parameter estimates and performance comparison (Gauging site: Punjnad).

DistributionMethodμ̂α̂k̂Total Rank
GLOLM7468.051 (671.275)2977.291 (317.703)-0.147 (0.07)67
MLE7431.049 (6.926)3144.326 (7.18)-0.285 (0.058)75
LS7341.052 (8.586)3466.268 (7.841)-0.224 (0.133)90
WLS7272.286 (5.094)3205.575 (4.879)-0.243 (0.103)86
RLS7164.378 (5.399)4421.462 (4.846)-0.6 (0.036)39
GEVLM5711.289 (602.06)4596.079 (467.391)0.035 (0.088)114
MLE5714.788 (8.1)4576.992 (8.463)0.026 (0.088)114
LS5719.329 (10.749)4604.797 (8.647)-0.031 (0.188)115
WLS5706.757 (6.21)4599.097 (5.75)-0.013 (0.139)119
RLS5699.688 (10.232)4591.782 (8.298)-0.235 (0.219)67
GUMLM5637.989 (581.505)4452.644 (487.011)---105
MLE5657.312 (11.528)4448.414 (10.605)---109
LS5625.473 (11.877)4462.63 (10.601)---110
WLS5659.205 (15.247)4466.814 (12.678)---114
RLS5641.932 (11.274)4451.137 (10.009)---106
GPDLM534.538 (424.809)11405.314 (2036.261)0.486 (0.144)143
MLE535.504 (5.419)11403.13 (5.824)0.48 (0.046)141
LS542.904 (10.416)11422.739 (9.598)0.49 (0.104)138
WLS533.236 (4.439)11405.137 (3.35)0.479 (0.07)143
RLS530.215 (11.619)11403.793 (8.518)0.444 (1.924)135
REVLM10778.26 (574.289)4452.644 (466.274)---27
MLE10774.825 (8.435)4446.246 (8.978)---28
LS10777.054 (11.633)4451.614 (10.348)---30
WLS10751.577 (14.371)4457.8 (11.693)---33
RLS10774.558 (9.962)4440.219 (12.089)---27
Table 9.

Estimated flood quantiles for different return years with 95% confidence intervals.

Non-exceedance probability (F) for different return years
Gauging Site: (Suited Distribution)Value0.80.90.960.980.990.9950.9980.999
5 years10 years25 years50 years100 years200 years500 years1000 years
Marala (GPD-LS)Upper17949.0427757.8245021.3462343.7384473.74112745.7162504.5212433
Fit14220.3719501.9726957.6933098.8539824.2547309.5158762.0669000.38
Lower11789.6714657.2317702.119548.1321075.222338.4323679.7724492.99
S.E1647.5873520.4467427.17911814.117913.4326421.1743318.9162844.6
Khanki (GPD-LS)Upper20667.4932509.6953739.8775428.17103573.4140097.8205577.2272469.6
Fit16128.1622377.3731285.0538689.146854.5655995.1570040.5582606.77
Lower13412.6516887.3420658.822992.3624956.5426609.8128400.7129508.81
S.E1919.5344116.6778694.00713784.120728.7730117.1147727.6266569.8
Qadirabad (GPD-WLS)Upper19765.3829280.2943891.5156717.9471304.9187894.08113368.6135731.4
Fit16104.7521707.1128789.433979.1539084.3244155.1350887.1956062.28
Lower13417.416626.2919779.421545.1422910.3323965.8525003.0125583.82
S.E1648.6713271.2526196.6869003.38312346.716269.4822453.127997.8
Trimmu (GEV-MLE)Upper13054.4118645.4928389.6638217.6150912.6667374.9696896.5127063
Fit12101.5116177.2222391.427947.2534442.3742088.254398.5765788.5
Lower11246.214105.8317799.6320564.2623344.6726151.7929931.232851.92
S.E462.7011165.0322735.4564580.7987199.15910847.8617859.1225457.34
Punjnad (GPD-WLS)Upper14811.3518987.1923236.5525703.9527670.5829238.0630832.9531758.61
Fit13463.0916678.7419634.6221184.722321.2723156.6423933.6924346.23
Lower12346.1214833.7416882.3517837.6418470.2418889.1419234.1819395.3
S.E640.88251080.5291653.5922047.2632394.8612694.4573020.9883221.168

[i] Fit: estimated flood quantiles;

[ii] S.E: bootstrap standard errors;

[iii] Lower (Upper): lower (upper) limits of 95% confidence intervals for estimated flood quantiles.

From these results, there is no single probability distribution that fits best for all gauging sites. This finding is in line with the other studies focused on at-site flood frequency analysis [5, 6, 22]. From the results, it is also clear that GP distribution with LS and WLS estimation is a suitable model for most of the sites. As there were varying sample sizes for different sites, it means that sample size does not have an impact on the performance of GP distribution. It is interesting to note that for all gauging sites, the most suitable probability model with one estimation method is also among the top performers with other estimation methods.

Once the best-fit probability distribution is chosen for any particular site, the next objective of the flood frequency analysis is to compute quantiles for different return periods. Such quantiles for different return periods are computed using quantile functions and parameter estimates of the most appropriate blend of probability distribution and method of estimation. For each gauging site, the quantile estimate xT with non-exceedance probability F for return periods 5, 10, 25, 50, 100, 200, 500, and 1000 years is given in Table 9. Uncertainties are always associated with the estimated flood quantiles computed for different return periods. To give an idea of these uncertainties, Table 9 also contains the 95% confidence intervals for estimated flood quantiles using parametric bootstrapping standard errors. For all gauging stations, the values of standard errors indicate that there is more uncertainty in flood estimates for longer return periods as compared to estimates of relatively short return periods.

Conclusions

The modeling of maximum annual discharge at different gauging sites of the Chenab River is done through different probability models blended with different estimation methods. To recommend the best-suited pair of probability models and methods of estimation for each gauging site, we compared the performance of GLO, GEV, GUM, GP, and REV. GUM distributions were estimated with the LM, ML, LS, WLS, and RLS methods. The comparison is based on the total rank of 7 different measures of fit. For each gauging site, the probability distribution with the highest total rank score is considered most appropriate.

Different probability distributions emerge as most suitable for flood frequency analysis at different sites. Generalized Pareto with LS estimation is the most suitable distribution for Marala and Khanki gauging sites. For gauging sites at Qadirabad and Punjnad, GP distribution with the WLS method is found to be the most suitable combination, while GEV distribution with the ML method performs better than all other distributions for the Trimmu gauging station. It can be concluded that GPD with LS and WLS estimation is the most suitable combination for most of the gauging stations (in fact, for all sites except for the gauging site at Trimmu). The results are well in line with many studies focused on at-site flood frequency analysis of rivers in different countries. The emergence of different best-suited probability distributions for different gauging sites has been reported in many other studies [5, 6, 9-11, 22]. Similarly, our finding about the larger uncertainty associated with longer return periods is common and reported in different studies [5, 6].

These results can be used to study and forecast flood levels, the management of water reservoirs, and the planning of hydraulic structures at the Chenab River (and related catchment areas). The most suitable distributions recommended for different gauging sites can be useful candidate distributions for regional analysis of the Chenab River or at-site analysis of other rivers in near geographical locations.

Limitations

Like every research project, the present study also has certain limitations, and work can be done in the future to overcome them. For instance, some other probability distributions, particularly four or five parameter distributions, can be applied. Similarly, different other estimation methods can be used for parameter estimation, like the method of maximum product spacing, Bayesian estimation, etc. Moreover, the analysis can be improved by incorporating a multivariate structure by taking temperature, rainfall, humidity, and other environmental factors into account.