The parameter estimation of statistical distributions is important for regional frequency analysis (RFA). The accuracy of different parts of RFA such as estimating the regional quantiles of the selected statistical distribution, determining the heterogeneity measure, and choosing the best distribution based on the Monte Carlo simulation, may be influenced by using the different values of regional parameters. To fulfill this aim, in the present study, a new model is developed for regional drought frequency analysis. This model utilizes the L-moments approach and the adjusted charged system search as an advanced meta-heuristic algorithm, in which some modifications on the equations of the algorithm are performed to improve its standard variant. The verification of the regional parameters estimated by the new methodology yields accurate results compared to other models. Furthermore, this study illustrates the usefulness of the robust discordancy measure against the classic one. For this purpose, different values of the subset factors (α) are utilized in the robust discordancy measure, and finally, the best value of subset factor is found equal to 0.8, which can accurately recognize discordant sites within the region.

## INTRODUCTION

Drought is a normal part of any climate (Keshavarz *et al.* 2013). Generally, droughts take place when a region receives consistently below average precipitation (Hossein *et al.* 2012). Compared with other natural hazards, the damaging effects of droughts are non-structural, and much larger in terms of the spatial extent (Reddy & Ganguli 2012), where 22% of economic damages and 33% of the number of affected persons are related to droughts (Belayneh *et al.* 2014). Therefore, utilizing regional frequency analysis (RFA) seems necessary to recognize the properties of this phenomenon. RFA can circumvent the limitation of statistical estimates, such as the absence of lengthy records or too short records (Zhang *et al.* 2012). RFA leads to more accurate quantile estimations than at-site analysis (Santillán *et al.* 2014).

To perform RFA, the considered region should be homogeneous, otherwise, it should be subdivided into some homogeneous regions (Seckin *et al.* 2013). Among many homogeneity tests in the hydrologic literature, the L-moments (LMOM)-based tests proposed by Hosking & Wallis (1993, 1997) are known as the most powerful ones (Viglione *et al.* 2007). Thus, numerous studies have utilized the LMOM approach to the RFA of hydrologic phenomena such as rainfall, stream flow, drought, and other fields of water engineering (Modarres 2010; Chérif & Bargaoui 2013; Dikbas *et al.* 2013; Rahman *et al.* 2013; Rajsekhar *et al.* 2013; Zakaria & Shabri 2013; Dodangeh *et al.* 2014; Feng *et al.* 2014; Goyal & Gupta 2014; Sarhadi & Heydarizadeh 2014).

Detection of the possible outlying data in a region is an important matter for RFA (Saf 2010). Hosking & Wallis (1993, 1997) recommended the discordancy measure based on the mean and covariance of the LMOM ratios using the Mahalanobis distance. This measure can be strongly affected by the discordant sites and is not robust for outlying data. To improve this problem, Neykov *et al.* (2007) proposed the robust discordancy measure. A few studies have then used this measure in hydrology. Saf (2010) analyzed the effect of discordancy detection measures on regional flood probability types and the accuracy of the estimation on regional analysis. Ilorme & Griffis (2013) utilized the robust discordancy measure to determine physically discordant sites and proposed a new methodology to identify the physical attributes that are the greatest demonstrator of extreme hydrological response. The mentioned studies showed the capability of robust discordancy measure against the classic one for outliers’ identification; however, the importance of subset size for all data sets was not investigated. In this study, the different subset sizes are considered to evaluate the accuracy of the robust discordancy measure.

The LMOM method is also utilized for parameter estimation of statistical distribution for RFA, although there are no simple expressions and explicit solutions for fulfilling this aim. Therefore, several approximations or simplifications are suggested for solving such nonlinear and complex problems (Hassanzadeh *et al.* 2011; Reddy & Singh 2014). Meta-heuristic algorithms (MHA) are one useful tool to solve these problems and various MHA have been extensively used to solve the optimization problems of hydrologic phenomena (e.g., Rai *et al.* 2009; Hassanzadeh *et al.* 2011; Reddy & Singh 2014; Shin *et al.* 2014).

Recently, an MNA known as charged system search (CSS) was introduced by Kaveh & Talatahari (2010), and applied to different engineering optimization problems. CSS utilizes the Coulomb and Gauss's laws from electrostatics, and the Newtonian laws of mechanics. There are few studies on applying CSS in the field of water resources. Kaveh *et al.* (2012) developed four different models of composite open channels and then optimized them by CSS. The results of this method were compared to those of ant colony optimization and genetic algorithm (GA). Talatahari *et al.* (2013) applied CSS and particle swarm optimization (PSO) for identifying the parameters of the linear Muskingum model. In order to evaluate these algorithms, a numerical example was utilized and the results were compared to those of other algorithms. Sheikholeslami *et al.* (2014) utilized CSS to optimal cost design of water distribution networks. They utilized some well-known benchmark instances as case studies. Comparison of the results of CSS with some other MHA demonstrated the high performance of this algorithm. In this study, an adjusted CSS (ACSS) algorithm is applied for estimating the regional parameters of Kappa distribution, and then the results are compared with those of the LMOM as the conventional approach.

Considering these points, the main contributions of this paper are:

- (1)
presenting the robust discordancy measure with different values of subset factor and comparing with the classic discordancy measure;

- (2)
applying the ACSS algorithm for the parameter estimation of statistical distributions; and

- (3)
proposing a new model to give regional parameters for statistical distribution.

## METHODOLOGY

### ACSS

The ACSS algorithm as a meta-heuristic algorithm contains a number of charged particles (CPs), where each one treats as a charged sphere and can insert an electric force to the others. Thus, the resultant electrical force affect a CP results in its acceleration. ACSS uses the governing laws of motion from the Newtonian mechanics to determine the position of CPs. Application of these laws provides a good balance between the exploration and the exploitation of the algorithm (Kaveh & Talatahari 2010). The pseudo-code for the ACSS algorithm is summarized as follows.

**Step 1: Initialization.**The magnitude of the charge for each CP is defined as: where

*f*and

_{b}*f*are the best and the worst objective function values among all of the particles;

_{w}*f*represents the fitness of the agent

_{i}*i*; and

*N*is the total number of CPs. The separation distance

*r*between any two CPs is defined as follows: where X

_{ij}*and X*

_{i}*are the positions of the*

_{j}*i*th and

*j*th CPs, respectively; X

*is the position of the best current CP; and*

_{b}*ɛ*is a small positive number. The initial positions of CPs are determined randomly and the initial velocities of CPs are assumed to be zero.

**Step 2: Charged memory creation.** A number of the best CPs and the values of their corresponding objective functions are saved in the charged memory.

**Step 3: Force determination.**The probability of moving each CP towards the others is determined using the following function: Then, the resultant force vector for each CP is calculated as: where

**F**

*is the resultant force acting on the*

_{j}*j*th CP; X

*and X*

_{i}*are the positions of the*

_{j}*i*th and

*j*th CPs, respectively.

**Step 4: Force type identification.**The kind of forces can be attractive or repelling, and this is determined by using the kind of force parameter,

*ar*, defined as: where

_{ij}*ar*determines the type of the force, with +1 representing the attractive force and −1 denoting the repelling force, and

_{ij}*k*is the force type coefficient, which controls the effect of the kind of force.

_{t}*ar*should be multiplied to the resultant force.

_{ij}**Step 5: Solution construction.**Each CP moves to the new position as: where

*k*and

_{a}*k*are the acceleration and the velocity coefficients, respectively;

_{v}*rand*and

_{j1}*rand*are two random numbers; V

_{j2}

_{j,}_{old}and V

_{j,}_{new}are the old and new velocities of the

*j*th CP, respectively; and Δ

*t*is the time step which is set to unity.

**Step 6: Velocities correction.**To prevent escaping CPs from the search area, the velocities of CPs in each dimension are restricted to a maximum velocity (v

_{j,}_{max}), as: The maximum velocity is defined as: where

*γ*is the fraction of the distance between the bounds in the range [0,1]. x

_{j}_{,max}and x

_{j}_{,min}represent the upper and lower bounds for each design variable, respectively.

**Step 7: Solutions modification.**The nearest limit and reflection methods is utilized for modifying the positions and velocities of violated CPs, respectively. Therefore, any variable of the violated CP is regenerated as follows: This step is utilized in the algorithm in order to direct the unfeasible solutions.

**Step 8: Updating of charged memory.** The best new vectors are included in the CM and the worst ones are excluded from the charged memory.

**Step 9: Terminating criterion control.** Steps 3–9 are repeated until a terminating criterion is satisfied.

It is notable that the steps 4, 6, and 7 are only applied in the ACSS algorithm.

### LMOM approach

Hosking (1986, 1990) defined LMOM as linear combinations of probability weighted moments, which can be interpreted as measures of the location, scale, and shape of probability distributions. Based on LMOM, Hosking & Wallis (1993) proposed three useful statistics in RFA, the discordancy, heterogeneity, and goodness-of-fit measures. Detailed explanations and formula can be found in Hosking & Wallis (1997).

### Discordancy measure

*D*) for the

_{i}*i*th site in a region is defined as: where

*N*is number of sites,

*u*= [

_{i}*t*

^{(i)},

*t*

_{3}

^{(i)},

*t*

_{4}

^{(i)}]

*is a vector of the LMOM ratios for the*

^{T}*i*th site, and its components are:

*t*as the L-coefficient of variation (L-CV),

*t*

_{3}as the L-coefficient of skewness (L-CS), and

*t*

_{4}as the L-coefficient of kurtosis (L-CK), respectively; is the regional unweighted average and

*S*is the matrix of sums of squares and cross-products. The critical discordancy measure is equal to 3 and the site becomes discordant when

*D*> 3 (Hosking & Wallis 1993, 1997). Due to the outlier observations’ influence on the classical estimates of and

_{i}*S*, the number of discordant sites for a region can be smaller than the real number. Therefore, it is necessary to utilize the robust discordancy measure in RFA (Saf 2010).

There are many robust estimators in statistics, such as minimum covariance determinant (MCD), minimum volume ellipsoid, and M-estimators (Alameddine *et al.* 2010). The MCD estimator is a highly robust estimator of multivariate location and scatter (Rousseeuw & Van Driessen 1999; Hubert *et al.* 2005; Rousseeuw *et al.* 2006). Rousseeuw (1984, 1985) defined the MCD estimator, which can be computed efficiently with the FAST-MCD algorithm. Finding the *h* observations (out of *N*) whose classical covariance matrix has the lowest possible determinant, is the main objective of MCD. There is a relation between the subset size (*h*) and subset factor (*α*) when the value of *α* varies from 0.5 to 1. If *α**=* 0.5, *h* equals , and when *α**=* 1, *h* equals *N* (Pison *et al.* 2002; Neykov *et al.* 2007). For other values of *α,* the *h* can be calculated by linear interpolation.

*et al.*(2007) introduced the robust discordancy measure (

*RD*) which uses the robust estimates (

_{i}*M*and

*C*), instead of classic ones ( and

*S*) in the Mahalanobis distance formula for RFA: where

*M*and

*C*are the average of

*h*observation and its covariance matrix (Pison

*et al.*2002; Hubert

*et al.*2008).

*w*is defined based on the initial

_{i}*RD*by: where

_{i}*χ*is a chi-square distribution with degrees of freedom by

^{2}*p*for the corresponding probability 0.975. For

*p*

*=*3, the critical discordancy value becomes The reweighted MCD estimators are calculated as: The final robust discordancy measure is computed using these reweighted estimates (

*M*,

_{R}*C*).

_{R}In this study, the MATLAB library for robust analysis, LIBRA, introduced by Verboven & Hubert (2005), is used for the computation of the robust discordancy measure.

### Heterogeneity measure

One of the main steps in the RFA contains identifying homogenous regions (Saf 2009). For this purpose, the heterogeneity measure was proposed to compute the degree of heterogeneity in a region. The steps of the heterogeneity measure are illustrated as below (Hosking & Wallis 1993).

**Step 1:**Calculate the regional LMOM ratios.**Step 2:**Compute the*V*variables.**Step 3:**Fit the four-parameter Kappa (KAP) distribution.**Step 4:**Simulate a large number of synthetic regions.

Hosking & Wallis (1993) suggested that a region is ‘acceptably homogeneous’ if *H _{i}* is less than 1, ‘possibly heterogeneous’ if

*H*is between 1 and 2, and ‘definitely heterogeneous’ if

_{i}*H*is greater than 2. Also, they stated that the

_{i}*H*statistic compared to the statistics

_{1}*H*

_{2}and

*H*

_{3}has a discriminating power and is the most important heterogeneity measure (Rao & Hamed 2000).

### Estimation of regional parameters

In order to determine the heterogeneity measure, it is necessary to estimate the regional parameters of KAP distribution for the sample data. In this study, three models are considered to fulfill this aim.

**Model 1:**In the first model, the parameters of the KAP distribution for each site are estimated using the LMOM and ACSS methods, separately. The cost function which is selected to be applied by the ACSS method in the

*i*th site is calculated as (Hassanzadeh

*et al.*2011): where is the

*j*th observed severity of the site

*i*; is the

*j*th computed severity of the site

*i*; and is the average of observed severity of the site

*i*.

**Model 2:** Based on this model, after combining the sites, the regional parameters of these sites are estimated using the proposed methods. For calculating the cost function, the combined data are used instead of at-site data.

**Model 3:**Finally, in the third model, all data of each site are used to estimate the regional parameters. In other words, the cost function is redefined in a way that the result will be matched with all sites with a small error, i.e.,

Due to the nature of the defined optimization problem, the regional parameters of KAP distribution can only be estimated by the ACSS method as a meta-heuristic algorithm.

### Goodness-of-fit measure

*Z*-based test (Hosking & Wallis 1993) and LMOM ratios-based test (Liou

^{DIST}*et al.*2008). The goodness-of-fit measure,

*Z*, which is utilized in this study for selecting the best regional distribution, is defined as: where is the L-CK of the fitted distribution to the data using the candidate distribution,

^{DIST}*DIST*can be any of three-parameter distributions, such as generalized Pareto (GPA), generalized extreme value (GEV), generalized logistic (GLO), lognormal (LN3) or Pearson type III (PE3). The value of for each probability distribution can be expressed by the following polynomial approximations (Wu

*et al.*2012):

The coefficient *A _{k}* was given by Hosking & Wallis (1997, Table A3).

*B*

_{4}and

*σ*

_{4}are the bias and the standard deviation of , respectively, as: where is the regional L-CK values for the

*m*th simulated region, and

*N*is the number of simulated regions with

_{sim}*N*sites. As mentioned in the heterogeneity measure, the KAP distribution is used in a similar way for the simulated region. For selected distribution, if , then it can be suitable as a regional distribution (Hosking & Wallis 1997).

### Estimation of the regional distribution

*N*sites, the quantile of nonexceedance probability

*F*for each site () is expressed as: where is the scale factor (mean of data at site

*i*), and is the regional growth curve. The dimensionless rescaled data , with

*i*

*=*1, …,

*N*and

*j*

*=*1, …,

*n*, are the basis for estimation . Then, the regional growth curve is computed using these parameters and the nonexceedance probability

_{i}*F*. Finally, Monte Carlo simulation is utilized to quantify the accuracy of the quantile estimates for various probabilities. Three different measures, such as the regional average relative bias (), regional average absolute relative bias (), and regional average relative RMSE () have been defined by Hosking & Wallis (1997) to fulfill this aim and here we utilized all of them in assessing the accuracy of quantile estimates by the following formulas: where is the site-

*i*quantile estimate for nonexceedance probability

*F*at the

*m*th simulated region.

### Cluster analysis

The cluster analysis is used to assemble objects into a set of specific groups with a maximum similarity between the cluster members (Modarres 2010). Generally, the clustering algorithms can be classified into two types: hierarchical and partitional (Demirel *et al.* 2012). Among the different partitional algorithms, the commonest and most well-known method is the K-means method (Sönmez & Kömüşcü 2011; Kar *et al.* 2012; Rahman *et al.* 2013).

*et al.*2011; Dikbas

*et al.*2012; Kar

*et al.*2012): where

*Y*is the normalized value, and

*X*and

_{min}*X*are the minimum and maximum values for each variable in the data set (

_{max}*X*), respectively.

*k*), and to check the success of cluster assigning, the cluster validity indices are extensively used in a data set (Rousseeuw 1987; Rao & Srinivas 2006). In the present study, the silhouette width is used for this aim as: where

*a*(

*i*) is the average dissimilarity of the

*i*th object with all other objects within the same cluster and

*b*(

*i*) is the lowest average dissimilarity of the

*i*th object with all other objects within any neighbor cluster. For the given

*k*clusters, the average silhouette width () is the average of the silhouette widths over all objects (Kaufman & Rousseeuw 1990).

### Adaption to droughts

Since droughts are regional in nature, it is logical to utilize RFA for drought assessment (Tallaksen *et al.* 1997). Various indices have been introduced to define drought characteristics, and commonly, the standardized precipitation index (SPI) is the most popular and widely used one for drought analysis (Raziei *et al.* 2009; Serinaldi *et al.* 2009; Santos *et al.* 2011; Mirabbasi *et al.* 2012; Chen *et al.* 2013). Among different time scales, the SPI values on a six-month time scale (SPI-6) is most useful for describing the shallow soil moisture available to agricultural crops (Reddy & Singh 2014; Abdi *et al.* 2016a, 2016b). After computing the SPI-6 values, the drought characteristics can be obtained from these values. In this study, we considered the drought severity (*S _{d}*) as one of the main characteristics of drought. The drought severity is defined based on the absolute value of cumulative SPI values within the duration of a drought event (Chen

*et al.*2013). The drought duration is defined as the number of consecutive intervals (months) where the SPI values are less than the threshold value of zero (Mirabbasi

*et al.*2012).

In order to adopt the LMOM approach to the drought severity, the details of the methodology can be summarized as follows:

Consider a number of sites within the region.

Assemble the monthly precipitation data of each site.

Calculate the SPI values and extract the drought characteristics.

Compute the discordancy and heterogeneity measures.

Specify the homogeneous sub-regions by using the cluster analysis.

Determine the best regional distribution for each sub-region.

### Study area and data

^{2}(Figure 1). Monthly precipitation data for a period of 31 years from 1982 to 2013 are extracted for 60 gauging sites from the East-Azarbaijan Meteorological Organization. The recorded minimum and maximum air temperature is −26 and +32 °C, respectively, with a weighted average of 10.2 °C (Zarghami

*et al.*2011). The climate of the region is classified as semi-arid cold climate based on the Amberje climate classification. The attributes at all sites in the study are shown in Table 1. The summary statistics of regional drought severity calculated for all sites based on the LMOM are (L-CV)

_{R}= 0.6548, (L-CS)

_{R}= 0.4807, and (L-CK)

_{R}= 0.2440.

Attribute | Statistic | Range |
---|---|---|

Annual precipitation | Average (mm) | 199.9 to 508.3 |

Maximum (mm) | 306.0 to 1,101.7 | |

Minimum (mm) | 70.1 to 287.3 | |

Coefficient of variation | 0.2 to 0.6 | |

Drought severity | Average | 4.2 to 6.4 |

Maximum | 14.5 to 84.7 | |

Minimum | 0.0 to 0.1 | |

Coefficient of variation | 0.7 to 2.9 | |

L-CV | 0.4 to 0.8 | |

L-CS | 0.1 to 0.8 | |

L-CK | −0.1 to 0.7 | |

Geographical coordination | Latitude (N) | 37.1 to 39.1 |

Longitude (E) | 45.4 to 47.9 |

Attribute | Statistic | Range |
---|---|---|

Annual precipitation | Average (mm) | 199.9 to 508.3 |

Maximum (mm) | 306.0 to 1,101.7 | |

Minimum (mm) | 70.1 to 287.3 | |

Coefficient of variation | 0.2 to 0.6 | |

Drought severity | Average | 4.2 to 6.4 |

Maximum | 14.5 to 84.7 | |

Minimum | 0.0 to 0.1 | |

Coefficient of variation | 0.7 to 2.9 | |

L-CV | 0.4 to 0.8 | |

L-CS | 0.1 to 0.8 | |

L-CK | −0.1 to 0.7 | |

Geographical coordination | Latitude (N) | 37.1 to 39.1 |

Longitude (E) | 45.4 to 47.9 |

## RESULTS

### Assumption of homogeneous region for all sites

First, the entire the study area is assumed as one homogeneous region. Hence, the classic discordancy (*CD _{i}*) measure and robust discordancy (

*RD*) measure with five different values of

_{i}*α*are computed and the discordant sites are shown in Table 2. As can be seen, sites 2, 43, 51, and 56 are the discordant sites according to the values of

*CD*and

_{i}*RD*. The difference between

_{i}*CD*and

_{i}*RD*with

_{i}*α*= 0.75, 0.95 is only at site 13, and with

*α*

*=*0.80, 0.85, 0.90, it is at sites 3 and 13.

Discordancy measure | |||||
---|---|---|---|---|---|

Robust | |||||

Classic | α = 0.75 | α = 0.80 | α = 0.85 | α = 0.90 | α = 0.95 |

2(3.45), 43(5.05), 51(3.24), 56(3.11)^{a} | 2(3.20), 13(3.67), 43(4.92), 51(3.88), 56(4.08) | 2(3.66), 3(3.06), 13(3.80), 43(4.83), 51(3.93), 56(3.85) | 2(3.82), 3(3.40), 13(4.10), 43(4.78), 51(4.18), 56(3.82) | 2(3.82), 3(3.40), 13(4.10), 43(4.78), 51(4.18), 56(3.82) | 2(3.53), 13(3.12), 43(4.87), 51(3.13), 56(3.73) |

Discordancy measure | |||||
---|---|---|---|---|---|

Robust | |||||

Classic | α = 0.75 | α = 0.80 | α = 0.85 | α = 0.90 | α = 0.95 |

2(3.45), 43(5.05), 51(3.24), 56(3.11)^{a} | 2(3.20), 13(3.67), 43(4.92), 51(3.88), 56(4.08) | 2(3.66), 3(3.06), 13(3.80), 43(4.83), 51(3.93), 56(3.85) | 2(3.82), 3(3.40), 13(4.10), 43(4.78), 51(4.18), 56(3.82) | 2(3.82), 3(3.40), 13(4.10), 43(4.78), 51(4.18), 56(3.82) | 2(3.53), 13(3.12), 43(4.87), 51(3.13), 56(3.73) |

^{a}Site number (*D _{i}* value).

Then, the regional parameters of the KAP distribution are estimated based on three proposed models. For model 1, the results of at-site parameters and the values of two main goodness-of-fit tests, the root mean square error (RMSE) and efficiency coefficient (EC), based on the LMOM and ACSS methods, are presented in Tables 3 and 4 for the selected sites, respectively. Tables 3 and 4 indicate that the results of the ACSS method are better than the LMOM method for all sites, especially for sites 13, 51, and 60 where the EC value becomes smaller than 0.9 for LMOM, whereas the performances of ACSS are adequate for these sites. Also, the values of *ξ* and *α* parameters of the LMOM results for sites 22, 42, 52, and 58 are significantly different in comparison with these parameters for the other sites. Therefore, the regional parameters can be extensively affected by the weighted ratio () of model 1. In contrast, the values of these parameters have some small variations for the results of the ACSS method.

Parameters of KAP distribution | Goodness-of-fit tests | |||||
---|---|---|---|---|---|---|

Site number | ξ | α | k | h | EC | RMSE |

13 | −0.3371 | 0.3914 | −0.6539 | 2.0483 | 0.8579 | 1.0119 |

22 | −15.8244 | 19.6832 | 0.9210 | 4.3378 | 0.9809 | 0.1882 |

42 | −21.7919 | 23.9107 | 0.8254 | 5.4338 | 0.9853 | 0.1879 |

51 | −0.3343 | 0.3368 | −0.6962 | 2.2837 | 0.8293 | 1.1716 |

52 | −41.4305 | 95.4353 | 2.1477 | 3.5142 | 0.9515 | 0.2266 |

58 | −16.3784 | 25.8720 | 1.2713 | 3.6308 | 0.9894 | 0.1179 |

60 | −3.9450 | 3.6644 | 0.2142 | 3.4678 | 0.8728 | 0.5099 |

Parameters of KAP distribution | Goodness-of-fit tests | |||||
---|---|---|---|---|---|---|

Site number | ξ | α | k | h | EC | RMSE |

13 | −0.3371 | 0.3914 | −0.6539 | 2.0483 | 0.8579 | 1.0119 |

22 | −15.8244 | 19.6832 | 0.9210 | 4.3378 | 0.9809 | 0.1882 |

42 | −21.7919 | 23.9107 | 0.8254 | 5.4338 | 0.9853 | 0.1879 |

51 | −0.3343 | 0.3368 | −0.6962 | 2.2837 | 0.8293 | 1.1716 |

52 | −41.4305 | 95.4353 | 2.1477 | 3.5142 | 0.9515 | 0.2266 |

58 | −16.3784 | 25.8720 | 1.2713 | 3.6308 | 0.9894 | 0.1179 |

60 | −3.9450 | 3.6644 | 0.2142 | 3.4678 | 0.8728 | 0.5099 |

Parameters of KAP distribution | Goodness-of-fit tests | |||||
---|---|---|---|---|---|---|

Site number | ξ | α | k | h | EC | RMSE |

13 | −0.3525 | 0.2005 | −1.0000 | 3.1951 | 0.9975 | 0.1334 |

22 | −4.7943 | 5.7087 | 0.5174 | 2.9539 | 0.9863 | 0.1595 |

42 | −7.4693 | 7.8411 | 0.5289 | 3.6934 | 0.9885 | 0.1658 |

51 | −0.7334 | 0.2095 | −1.0000 | 4.9999 | 0.9988 | 0.0996 |

52 | −3.9892 | 7.5960 | 1.0636 | 2.0900 | 0.9585 | 0.2096 |

58 | −5.9095 | 8.9804 | 0.9055 | 2.6776 | 0.9897 | 0.1163 |

60 | −8.9238 | 9.9841 | 0.6695 | 3.9585 | 0.9314 | 0.3744 |

Parameters of KAP distribution | Goodness-of-fit tests | |||||
---|---|---|---|---|---|---|

Site number | ξ | α | k | h | EC | RMSE |

13 | −0.3525 | 0.2005 | −1.0000 | 3.1951 | 0.9975 | 0.1334 |

22 | −4.7943 | 5.7087 | 0.5174 | 2.9539 | 0.9863 | 0.1595 |

42 | −7.4693 | 7.8411 | 0.5289 | 3.6934 | 0.9885 | 0.1658 |

51 | −0.7334 | 0.2095 | −1.0000 | 4.9999 | 0.9988 | 0.0996 |

52 | −3.9892 | 7.5960 | 1.0636 | 2.0900 | 0.9585 | 0.2096 |

58 | −5.9095 | 8.9804 | 0.9055 | 2.6776 | 0.9897 | 0.1163 |

60 | −8.9238 | 9.9841 | 0.6695 | 3.9585 | 0.9314 | 0.3744 |

The regional parameters of KAP distribution based on model 1 are given in Table 5. The results of the goodness-of-fit tests by using the regional parameters of model 1 for each site demonstrate that LMOM performs worse than ACSS. For 97% of sites, the results of EC are less than zero for LMOM while for ACSS only 8% of sites have EC < 0 and in return 65% and 30% of sites have EC greater than 0.5 and 0.8, respectively.

Parameters of the KAP distribution | |||||
---|---|---|---|---|---|

Model number | Method | ξ | α | K | h |

1 | LMOM | −3.3795 | 5.0803 | 0.1223 | 2.3682 |

ACSS | −2.3047 | 2.9404 | 0.0998 | 2.338 | |

2 | LMOM | −1.0803 | 1.4655 | −0.0731 | 2.0612 |

ACSS | −1.0179 | 1.4671 | −0.0613 | 1.9648 | |

3 | ACSS | −0.7107 | 1.4442 | 0.0256 | 1.6082 |

CSS | −0.6209 | 1.3793 | 0.0109 | 1.5196 | |

PSO | −0.4925 | 1.2782 | −0.0167 | 1.3974 | |

GA | −0.5197 | 1.2441 | −0.0371 | 1.4939 |

Parameters of the KAP distribution | |||||
---|---|---|---|---|---|

Model number | Method | ξ | α | K | h |

1 | LMOM | −3.3795 | 5.0803 | 0.1223 | 2.3682 |

ACSS | −2.3047 | 2.9404 | 0.0998 | 2.338 | |

2 | LMOM | −1.0803 | 1.4655 | −0.0731 | 2.0612 |

ACSS | −1.0179 | 1.4671 | −0.0613 | 1.9648 | |

3 | ACSS | −0.7107 | 1.4442 | 0.0256 | 1.6082 |

CSS | −0.6209 | 1.3793 | 0.0109 | 1.5196 | |

PSO | −0.4925 | 1.2782 | −0.0167 | 1.3974 | |

GA | −0.5197 | 1.2441 | −0.0371 | 1.4939 |

The regional parameters of KAP distribution for model 2 are given in Table 5. The results of the goodness-of-fit tests obtained by the regional parameters show that the performance of the LMOM and ACSS methods is close. However, the value of the sum of differences of the EC (SDEC) measure is equal to 0.4631, which shows the relative superiority of ACSS. In both methods, the value of EC for site 43 is less than zero. As shown in Table 2, site 43 is a discordant site and has the maximum value of the discordancy measure among the others. Therefore, this model is less sensitive than the first one to the outliers. In addition, the results demonstrate that the percentage of sites which have EC values greater than 0.5 and 0.8, obtained by LMOM are equal to 97% and 72%, respectively. These values are 97% and 75% for ACSS, respectively.

Table 5 presents the regional parameters of the KAP distribution estimated from model 3. Based on these parameters, the value of the SDEC measure for this model is equal to 1.8301. The results of the ACSS method for this model are more accurate than those for the second model. The value of EC in model 3 for site 43 is equal to 0.3345, whereas less than zero was obtained for the previous model. Indeed, this indicates that the sensitivity of model 3 to the discordant sites is less than model 2. The values of EC obtained by the ACSS results for model 3 show that 98% and 87% of the sites have values greater than 0.5 and 0.8, respectively. Consequently, model 3 yields the global optimal parameters, while the optimum parameters of model 2 are local.

Method | Best value of cost function | Mean of EC | Mean of RMSE |
---|---|---|---|

ACSS | 7.6924 | 0.8718 | 0.5023 |

CSS | 7.6958 | 0.8717 | 0.5028 |

PSO | 7.7079 | 0.8715 | 0.5029 |

GA | 7.7391 | 0.8710 | 0.5033 |

Method | Best value of cost function | Mean of EC | Mean of RMSE |
---|---|---|---|

ACSS | 7.6924 | 0.8718 | 0.5023 |

CSS | 7.6958 | 0.8717 | 0.5028 |

PSO | 7.7079 | 0.8715 | 0.5029 |

GA | 7.7391 | 0.8710 | 0.5033 |

According to the results of various models and methods, model 3 with the ACSS method is selected for estimating the regional parameters of the KAP distribution. However, with regional parameters of the KAP distribution, sample sizes the same as at-site historical data and *N _{sim}* = 1,000, the heterogeneity and goodness-of-fit measures are calculated as shown in Table 7, which indicate the study area with 60 sites is ‘definitely heterogeneous’ before removal of the discordant sites according to the

*H*

_{1}value. Also, based on the results of goodness-of-fit measure (

*Z*), only the PE3 distribution is selected for the study area.

^{DIST}After removing discordant sites based on discordancy measure | ||||||||
---|---|---|---|---|---|---|---|---|

Measure | Before removing discordant sites | Classic | Robust | |||||

α = 0.75 | α = 0.80 | α = 0.85 | α = 0.90 | α = 0.95 | ||||

Heterogeneity | H_{1} | 3.3 | 2.13 | 1.51 | 1.29 | 1.29 | 1.29 | 1.51 |

H_{2} | 5.06 | 4.12 | 3.49 | 3.33 | 3.33 | 3.33 | 3.49 | |

H_{3} | 6.68 | 5.58 | 4.66 | 4.53 | 4.53 | 4.53 | 4.66 | |

Goodness-of-fit | GPA | −3.8 | −4.71 | −4.68 | −4.62 | −4.62 | −4.62 | −4.68 |

GEV | −7.28 | −8.54 | −8.48 | −8.36 | −8.36 | −8.36 | −8.48 | |

GLO | −8.01 | −9.73 | −9.3 | −9.18 | −9.18 | −9.18 | −9.3 | |

LN3 | −4.39 | −5.64 | −5.4 | −5.38 | −5.38 | −5.38 | −5.4 | |

PE3 | 0.51 | 0.53 | −0.45 | −0.31 | −0.31 | −0.31 | −0.45 |

After removing discordant sites based on discordancy measure | ||||||||
---|---|---|---|---|---|---|---|---|

Measure | Before removing discordant sites | Classic | Robust | |||||

α = 0.75 | α = 0.80 | α = 0.85 | α = 0.90 | α = 0.95 | ||||

Heterogeneity | H_{1} | 3.3 | 2.13 | 1.51 | 1.29 | 1.29 | 1.29 | 1.51 |

H_{2} | 5.06 | 4.12 | 3.49 | 3.33 | 3.33 | 3.33 | 3.49 | |

H_{3} | 6.68 | 5.58 | 4.66 | 4.53 | 4.53 | 4.53 | 4.66 | |

Goodness-of-fit | GPA | −3.8 | −4.71 | −4.68 | −4.62 | −4.62 | −4.62 | −4.68 |

GEV | −7.28 | −8.54 | −8.48 | −8.36 | −8.36 | −8.36 | −8.48 | |

GLO | −8.01 | −9.73 | −9.3 | −9.18 | −9.18 | −9.18 | −9.3 | |

LN3 | −4.39 | −5.64 | −5.4 | −5.38 | −5.38 | −5.38 | −5.4 | |

PE3 | 0.51 | 0.53 | −0.45 | −0.31 | −0.31 | −0.31 | −0.45 |

*α*confirm this concept, whereas according to the classic discordancy measure, only site 51 is discordant. Consequently, it seems that the obtained number of discordant sites according to the classic discordancy measure is not exactly correct.

At last, after removing the discordant sites based on the classic and robust discordancy measures, the heterogeneity and goodness-of-fit measures are computed as demonstrated in Table 7. A value of *H*_{1} = 2.13 in this table indicates that the study area is ‘definitely heterogeneous’. On the contrary, using the robust discordancy measure with different values of subset factor for recognizing and removing the discordant sites shows an improvement in heterogeneity measure and suggests that the whole region is ‘possibly heterogeneous’. Also, based on the goodness-of-fit measure, the PE3 distribution is selected as the regional distribution for this region.

### Clustering and identification of homogeneous regions

In order to identify the optimum number of clusters, the average silhouette width for some number of clusters are computed and illustrated in Table 8. According to this table, the maximum value of average silhouette width corresponds to the best number of clusters (*k*), which is equal to 3.

Number of clusters (k) | Average silhouette width |
---|---|

2 | 0.4315 |

3 | 0.5122 |

4 | 0.4907 |

5 | 0.4704 |

6 | 0.4815 |

7 | 0.5050 |

Number of clusters (k) | Average silhouette width |
---|---|

2 | 0.4315 |

3 | 0.5122 |

4 | 0.4907 |

5 | 0.4704 |

6 | 0.4815 |

7 | 0.5050 |

The values in bold indicate the optimum number of cluster and its average silhouette width.

*k*= 3 are presented in Figure 5. The silhouette index for a few sites in cluster 2 is less than zero, so the quality of clustering may not be good. Thus, the exact results can be achieved by using the heterogeneity measure for each sub-region. Figure 6 shows the location of the sites in the identified clusters.

The classic discordancy measure and robust discordancy measure with five different values of *α* are computed and the discordant sites for three clusters are presented in Table 9. The heterogeneity and goodness-of-fit measures are demonstrated in Table 10, as well. The results in Table 10 show that the third cluster is ‘acceptably homogeneous’, the first cluster is ‘possibly heterogeneous’, and the second one is ‘definitely heterogeneous’ according to the *H*_{1} test.

Discordancy measure | ||||||
---|---|---|---|---|---|---|

Robust | ||||||

Cluster | Classic | α = 0.75 | α = 0.80 | α = 0.85 | α = 0.90 | α = 0.95 |

1 | – | 23, 37, 43 | 37, 43 | 37, 43 | 43 | 43 |

2 | 56 | 13, 51, 56 | 13, 51, 53, 56 | 13, 51, 56 | 13, 51, 56 | 53, 56 |

3 | – | 42, 48 | 42, 48 | 42, 48 | 42, 48 | 42 |

Discordancy measure | ||||||
---|---|---|---|---|---|---|

Robust | ||||||

Cluster | Classic | α = 0.75 | α = 0.80 | α = 0.85 | α = 0.90 | α = 0.95 |

1 | – | 23, 37, 43 | 37, 43 | 37, 43 | 43 | 43 |

2 | 56 | 13, 51, 56 | 13, 51, 53, 56 | 13, 51, 56 | 13, 51, 56 | 53, 56 |

3 | – | 42, 48 | 42, 48 | 42, 48 | 42, 48 | 42 |

After removing discordant sites based on discordancy measure | |||||||||
---|---|---|---|---|---|---|---|---|---|

Classic | Robust | ||||||||

Measure | Cluster | Before removing discordant sites | α = 0.75 | α = 0.80 | α = 0.85 | α = 0.90 | α = 0.95 | ||

Heterogeneity | H_{1} | 1 | 1.58 | 1.58 | −0.93 | −0.70 | −0.70 | 0.29 | 0.29 |

2 | 2.14 | 2.03 | 1.37 | 0.87 | 1.37 | 1.37 | 1.63 | ||

3 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | ||

H_{2} | 1 | 2.09 | 2.09 | 0.50 | 0.59 | 0.59 | 1.24 | 1.24 | |

2 | 3.44 | 3.55 | 2.75 | 2.19 | 2.75 | 2.75 | 3.13 | ||

3 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | ||

H_{3} | 1 | 3.25 | 3.25 | 1.90 | 1.72 | 1.72 | 2.36 | 2.36 | |

2 | 5.10 | 5.31 | 3.93 | 3.49 | 3.93 | 3.93 | 5.09 | ||

3 | 1.10 | 1.10 | 1.10 | 1.10 | 1.10 | 1.10 | 1.10 | ||

Goodness-of-fit | GPA | 1 | −2.48 | −2.48 | −2.52 | −2.58 | −2.58 | −2.52 | −2.52 |

2 | −2.88 | −3.25 | −3.66 | −3.71 | −3.66 | −3.66 | −3.13 | ||

3 | −1.93 | −1.93 | −1.93 | −1.93 | −1.93 | −1.93 | −1.93 | ||

GEV | 1 | −5.15 | −5.15 | v4.81 | −4.86 | −4.86 | −4.81 | −4.81 | |

2 | −5.07 | −5.50 | −6.00 | −5.96 | −6.00 | −6.00 | −5.17 | ||

3 | −3.84 | −3.84 | −3.84 | −3.84 | −3.84 | −3.84 | −3.84 | ||

GLO | 1 | −5.88 | −5.88 | −5.46 | −5.49 | −5.49 | −5.40 | −5.40 | |

2 | −5.44 | −5.88 | −6.44 | −6.37 | −6.44 | −6.44 | −5.49 | ||

3 | −4.26 | −4.26 | −4.26 | −4.26 | −4.26 | −4.26 | −4.26 | ||

LN3 | 1 | −3.57 | −3.57 | −3.50 | −3.50 | −3.50 | −3.33 | −3.33 | |

2 | −2.87 | −3.23 | −3.83 | −3.79 | −3.83 | −3.83 | −3.02 | ||

3 | −2.33 | −2.33 | −2.33 | −2.33 | −2.33 | −2.33 | −2.33 | ||

PE3 | 1 | −0.87 | −0.87 | −1.25 | −1.19 | −1.19 | −0.80 | −0.80 | |

2 | 0.83 | 0.61 | −0.15 | −0.13 | −0.15 | −0.15 | 0.60 | ||

3 | 0.23 | 0.23 | 0.23 | 0.23 | 0.23 | 0.23 | 0.23 |

After removing discordant sites based on discordancy measure | |||||||||
---|---|---|---|---|---|---|---|---|---|

Classic | Robust | ||||||||

Measure | Cluster | Before removing discordant sites | α = 0.75 | α = 0.80 | α = 0.85 | α = 0.90 | α = 0.95 | ||

Heterogeneity | H_{1} | 1 | 1.58 | 1.58 | −0.93 | −0.70 | −0.70 | 0.29 | 0.29 |

2 | 2.14 | 2.03 | 1.37 | 0.87 | 1.37 | 1.37 | 1.63 | ||

3 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | ||

H_{2} | 1 | 2.09 | 2.09 | 0.50 | 0.59 | 0.59 | 1.24 | 1.24 | |

2 | 3.44 | 3.55 | 2.75 | 2.19 | 2.75 | 2.75 | 3.13 | ||

3 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | ||

H_{3} | 1 | 3.25 | 3.25 | 1.90 | 1.72 | 1.72 | 2.36 | 2.36 | |

2 | 5.10 | 5.31 | 3.93 | 3.49 | 3.93 | 3.93 | 5.09 | ||

3 | 1.10 | 1.10 | 1.10 | 1.10 | 1.10 | 1.10 | 1.10 | ||

Goodness-of-fit | GPA | 1 | −2.48 | −2.48 | −2.52 | −2.58 | −2.58 | −2.52 | −2.52 |

2 | −2.88 | −3.25 | −3.66 | −3.71 | −3.66 | −3.66 | −3.13 | ||

3 | −1.93 | −1.93 | −1.93 | −1.93 | −1.93 | −1.93 | −1.93 | ||

GEV | 1 | −5.15 | −5.15 | v4.81 | −4.86 | −4.86 | −4.81 | −4.81 | |

2 | −5.07 | −5.50 | −6.00 | −5.96 | −6.00 | −6.00 | −5.17 | ||

3 | −3.84 | −3.84 | −3.84 | −3.84 | −3.84 | −3.84 | −3.84 | ||

GLO | 1 | −5.88 | −5.88 | −5.46 | −5.49 | −5.49 | −5.40 | −5.40 | |

2 | −5.44 | −5.88 | −6.44 | −6.37 | −6.44 | −6.44 | −5.49 | ||

3 | −4.26 | −4.26 | −4.26 | −4.26 | −4.26 | −4.26 | −4.26 | ||

LN3 | 1 | −3.57 | −3.57 | −3.50 | −3.50 | −3.50 | −3.33 | −3.33 | |

2 | −2.87 | −3.23 | −3.83 | −3.79 | −3.83 | −3.83 | −3.02 | ||

3 | −2.33 | −2.33 | −2.33 | −2.33 | −2.33 | −2.33 | −2.33 | ||

PE3 | 1 | −0.87 | −0.87 | −1.25 | −1.19 | −1.19 | −0.80 | −0.80 | |

2 | 0.83 | 0.61 | −0.15 | −0.13 | −0.15 | −0.15 | 0.60 | ||

3 | 0.23 | 0.23 | 0.23 | 0.23 | 0.23 | 0.23 | 0.23 |

In order to improve the heterogeneity measure, Table 9 demonstrates that the classic discordancy measure could not detect any discordant site for cluster 1. Also, only one discordant site is detected for cluster 2, whereas the different discordant sites are identified based on the robust discordancy measure with different values of *α*. Finally, after removing the discordant sites for clusters 1 and 2, the heterogeneity and goodness-of-fit measures are also shown in Table 10. In this condition, the value of *H*_{1} demonstrates that cluster 2 is not ‘acceptably homogeneous’. On the contrary, after removing the discordant sites (based on robust discordancy measure with all values of *α*), cluster 1 is found to be ‘acceptably homogeneous’. Just *α* = 0.80 yields ‘acceptably homogeneous’ for cluster 2 and the other values of *α* yield ‘possibly heterogeneous’.

After the selection of PE3 distribution as a suitable regional frequency distribution, the regional parameters of the respective distribution are estimated using model 3 and by the ACSS method, then the regional quantiles of various return periods for all clusters are illustrated in Table 11.

Cluster | Regional parameters | Regional quantiles | ||||||
---|---|---|---|---|---|---|---|---|

γ | α | β | T = 10 | T = 25 | T = 50 | T = 100 | T = 200 | |

1 | 1.0000 | 2.7212 | 0.1482 | 9.92 | 15.88 | 21.52 | 28.01 | 35.26 |

2 | 1.0000 | 3.2201 | 0.1524 | 11.61 | 19.12 | 26.21 | 34.33 | 43.38 |

3 | 1.0000 | 3.1939 | 0.1461 | 12.04 | 19.86 | 27.28 | 35.84 | 45.41 |

Cluster | Regional parameters | Regional quantiles | ||||||
---|---|---|---|---|---|---|---|---|

γ | α | β | T = 10 | T = 25 | T = 50 | T = 100 | T = 200 | |

1 | 1.0000 | 2.7212 | 0.1482 | 9.92 | 15.88 | 21.52 | 28.01 | 35.26 |

2 | 1.0000 | 3.2201 | 0.1524 | 11.61 | 19.12 | 26.21 | 34.33 | 43.38 |

3 | 1.0000 | 3.1939 | 0.1461 | 12.04 | 19.86 | 27.28 | 35.84 | 45.41 |

The estimated regional growth curve values and the relative RMSE for each cluster are determined in Table 12. The results verify that by increasing the value of the return period, the accuracy of the regional growth curve decreases. In addition, Monte Carlo simulation is used to assess the accuracy of the estimates of the determined homogeneous regions, by generating the same sample sizes as those of the original sample. The number of random samples, *N _{sim}*, is set to 10,000 for selected quantiles corresponding to return periods of 10, 25, 50, 100, and 200 years. Table 13 shows the results of the

*B*(

^{R}*F*),

*A*(

^{R}*F*), and

*R*(

^{R}*F*). According to Tables 12 and 13, the values of the RMSE of the estimated growth curve values are always less than the estimated quantiles for the same return periods, which represent the uncertainty in the estimated scale factor. These results are verified by Saf (2009) and Hussain (2011).

Cluster | Regional growth curve | RMSE | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

T = 10 | T = 25 | T = 50 | T = 100 | T = 200 | T = 10 | T = 25 | T = 50 | T = 100 | T = 200 | |

1 | 2.01 | 3.22 | 4.37 | 5.69 | 7.16 | 0.1174 | 0.1582 | 0.1820 | 0.2021 | 0.2196 |

2 | 2.25 | 3.70 | 5.07 | 6.65 | 8.40 | 0.1464 | 0.1951 | 0.2252 | 0.2516 | 0.2752 |

3 | 2.17 | 3.57 | 4.91 | 6.45 | 8.17 | 0.1368 | 0.1677 | 0.1792 | 0.1860 | 0.1901 |

Cluster | Regional growth curve | RMSE | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

T = 10 | T = 25 | T = 50 | T = 100 | T = 200 | T = 10 | T = 25 | T = 50 | T = 100 | T = 200 | |

1 | 2.01 | 3.22 | 4.37 | 5.69 | 7.16 | 0.1174 | 0.1582 | 0.1820 | 0.2021 | 0.2196 |

2 | 2.25 | 3.70 | 5.07 | 6.65 | 8.40 | 0.1464 | 0.1951 | 0.2252 | 0.2516 | 0.2752 |

3 | 2.17 | 3.57 | 4.91 | 6.45 | 8.17 | 0.1368 | 0.1677 | 0.1792 | 0.1860 | 0.1901 |

Measure | Cluster | Return period (years) | ||||
---|---|---|---|---|---|---|

T = 10 | T = 25 | T = 50 | T = 100 | T = 200 | ||

B(^{R}F) | 1 | −0.1280 | −0.2256 | −0.2604 | −0.2791 | −0.2892 |

2 | −0.1583 | −0.2580 | −0.2906 | −0.3070 | −0.3150 | |

3 | −0.1453 | −0.2601 | −0.3026 | −0.3271 | −0.3420 | |

A(^{R}F) | 1 | 0.1435 | 0.2280 | 0.2604 | 0.2791 | 0.2892 |

2 | 0.1634 | 0.2629 | 0.3014 | 0.3243 | 0.3408 | |

3 | 0.1540 | 0.2601 | 0.3026 | 0.3271 | 0.3420 | |

R(^{R}F) | 1 | 0.2056 | 0.2753 | 0.3052 | 0.3237 | 0.3358 |

2 | 0.2306 | 0.3059 | 0.3384 | 0.3591 | 0.3730 | |

3 | 0.2311 | 0.3067 | 0.3405 | 0.3611 | 0.3739 |

Measure | Cluster | Return period (years) | ||||
---|---|---|---|---|---|---|

T = 10 | T = 25 | T = 50 | T = 100 | T = 200 | ||

B(^{R}F) | 1 | −0.1280 | −0.2256 | −0.2604 | −0.2791 | −0.2892 |

2 | −0.1583 | −0.2580 | −0.2906 | −0.3070 | −0.3150 | |

3 | −0.1453 | −0.2601 | −0.3026 | −0.3271 | −0.3420 | |

A(^{R}F) | 1 | 0.1435 | 0.2280 | 0.2604 | 0.2791 | 0.2892 |

2 | 0.1634 | 0.2629 | 0.3014 | 0.3243 | 0.3408 | |

3 | 0.1540 | 0.2601 | 0.3026 | 0.3271 | 0.3420 | |

R(^{R}F) | 1 | 0.2056 | 0.2753 | 0.3052 | 0.3237 | 0.3358 |

2 | 0.2306 | 0.3059 | 0.3384 | 0.3591 | 0.3730 | |

3 | 0.2311 | 0.3067 | 0.3405 | 0.3611 | 0.3739 |

## CONCLUSIONS

In the standard CSS algorithm, search agents or CPs affect each other according to laws from electrostatics and mechanics. Here, a new ACSS is developed to improve the performance of CSS. The ACSS and CSS algorithms use the governing laws from electrical physics and Newtonian mechanics to determine the amount and the direction of a CPs movement. The CSS and ACSS algorithms can distinguish the local search and the global search phases, and utilize suitable relationships in these phases resulting in a good balance between exploration and exploitation. However, applying some controlling limits and the constraints is the main difference between the ACSS and CSS algorithms.

On the other hand, a new model is introduced to estimate the regional parameters of the statistical distributions. Then, the ACSS method is utilized to estimate the accurate parameters and quantile estimates for a region. Also for the Monte Carlo simulation, this new model is employed to generate samples from the selected distribution. In addition, the robust discordancy measure with different subset factors is developed to detect and remove the discordant sites, which causes some improvements in the heterogeneity measure. The data used in this study are drought severity extracted for 60 sites in East-Azarbaijan, Iran. The relevant conclusions can be summarized as follows:

The ACSS has the ability to improve the search and convergence speed compared to the standard CSS.

The ACSS algorithm has stable convergence characteristics and good computational ability, and is an effective algorithm compared to other algorithms, such as the GA and PSO, for RFA.

The ACSS can estimate the global optimum (or at least near it) of parameters without using nonlinear complex equations and any approximations.

The new model is the best model among the others. This model is based on the participation of all the data of sites in a region for estimating the regional parameters and quantiles, while other models are based on the weighted at-site parameters or the combination of the sites.

Identifying discordant sites using the robust discordancy measure is more accurate than the classic discordancy measure. Among the different values for the subset factor, the value of 0.8 is found suitable for robust discordancy measure and to detect the discordant sites of a heterogeneous region.

The study area is divided into three sub-regions, south, northeast, and northwest by using the K-means clustering method and the Pearson type 3 distribution is found the best regional distribution for all sub-regions.

Finally, this method can be applied for RFA of other hydrological phenomena, such as rainfall and stream flow. Also, the multivariate RFA framework based on the multivariate LMOM approach is suggested for future application of this work. In addition, another goodness-of-fit measure, the LMOM ratios-based test, can be applied and compared with the *Z ^{DIST}*-based test for selecting the best regional distribution.