Application of Harmony Search to Design Storm Estimation from Probability Distribution Models

Abstract

The precision of design storm estimation depends on the selection of an appropriate probability distribution model (PDM) and parameter estimation techniques. Generally, estimated parameters for PDMs are provided based on the method of moments, probability weighted moments, and maximum likelihood (ML). The results using ML are more reliable than the other methods. However, the ML is more laborious than the other methods because an iterative numerical solution must be used. In the meantime, metaheuristic approaches have been developed to solve various engineering problems. A number of studies focus on using metaheuristic approaches for estimation of hydrometeorological variables. Applied metaheuristic approaches offer reliable solutions but use more computation time than derivative-based methods. Therefore, the purpose of the current study is to enhance parameter estimation of PDMs for design storms using a recently developed metaheuristic approach known as a harmony search (HS). The HS is compared to the genetic algorithm (GA) and ML via simulation and case study. The results of this study suggested that the performance of the GA and HS was similar and showed more accurate results than that of the ML. Furthermore, the HS required less computation time than the GA.

1. Introduction

Floods that occurred due to severe rain or major storm generally result in damage to properties and negative impacts on human activity.  Flood estimation is the process of minimizing property damage and reducing the threat to human activity. Design storm based on precipitation frequency analysis is one of the main procesess for flood estimation as well as a statistical representation of a precipitation event.

The primary purpose of precipitation frequency analysis in hydrometeorology is to estimate the magnitude of a storm event with a given frequency of occurrence. It can also be used to estimate the frequency of occurrence of a storm event with a given magnitude [1]. The precision of precipitation frequency analysis depends on the selection of an appropriate probability distribution model (PDM) and parameter estimation techniques. A number of PDMs have been developed to describe the probability distribution of the hydrometeorological variables. In practice, it is often assumed that the correct PDM is a member of the developed PDMs. However, the selection of an appropriate PDM is still one of the major problems in precipitation frequency analysis [2].

For each of the developed PDMs, estimated parameters are provided based on alternative estimation techniques, such as the method of moments (MOM), probability weighted moments (PWM), linear function of ranked observations (L-moments), and maximum likelihood (ML) [1–5]. The MOM is one of the most simple and commonly used methods for estimating the parameters. PWM and L-moments are discussed by Stedinger et al. [6]. Estimates for PWM can be obtained from order statics. Additionally, Stedinger et al. [6] recommended a parameter estimator for PWM in regionalization studies. L-moments tend to produce less variable estimates for higher moments, when an unusually large or small observation happens to be present in a sample [1]. The moment estimators, including MOM, PWM, and L-moments, are a simpler way to obtain the PDM parameters but are less accurate than the ML estimators. The results using ML are generally more reliable than the other methods. However, the ML is more laborious than the other methods because an iterative numerical solution, such as the Newton-Raphson method [7], must be used for estimating the parameters.

Metaheuristic approaches have been developed to solve various engineering optimization problems (e.g., linear and stochastic, dynamic, and nonlinear). These approaches include genetic algorithms [8, 9], ant colony optimization [10], simulated annealing [11], tabu searches [12], and evolutionary computation methods. Metaheuristic approaches use a stochastic random search instead of a gradient search so that intricate derivative information is unnecessary [13]. Therefore, the metaheuristic approaches have been shown to be a useful strategy to solve optimization problems in hydrology [5, 14–21].

Karahan et al. [22] used the genetic algorithm (GA) to predict rainfall intensities for a given set of return periods. The results showed that the proposed GA could be used to develop rainfall intensity-duration-frequency relationships with the lowest mean-squared error between the observed and predicted intensities. Karahan et al. [22] concluded that predicted intensities were in good agreement with the analyzed return period. Hassanzadeh et al. [23] studied the most suitable PDM for annual maximum discharge in East-Azerbaijan, Iran. Hassanzadeh et al. [23] employed the GA and ant colony optimization (ACO) techniques to find the parameters of PDMs. The performance of these algorithms was evaluated by comparison with conventional methods, such as MOM, ML, and PWM. The results showed that the GA and ACO were effective optimization tools compared to other methods for the parameter estimation of PDMs. The GA and ACO techniques could also be used for systems that are more complex and involve nonlinear optimization problems. However, the GA and ACO techniques need more computation time than the MOM, ML, and PWM methods to find the parameters of PDMs. Although the GA and ACO techniques are reasonable alternatives to solve hydrological problems, large computation time is an obvious disadvantage.

Therefore, the purpose of this study is to enhance the design storm estimation with improving parameter estimation of PDMs using a recently developed metaheuristic approach, a harmony search (HS) by Geem et al. [24]. The performance of the HS approach is compared with the GA and conventional methods (i.e., ML).

This paper is organized in the following manner. Section 2 introduces the methodology for the parameters of PDMs with statistical test criteria. The results of the simulation are shown in Section 3, and a case study is reported in Section 4. Finally, the summary and conclusions are presented in Section 5.

2. Methodology

2.1. Probability Distribution Models

To test the performance of the proposed parameter estimation approach, the two-parameter lognormal (LN2), two-parameter gamma (GAM2), generalized extreme value (GEV), and Gumbel (GUM) distribution models are used. The general properties of the PDMs are given in Appendices.

2.2. Metaheuristic Approaches

In the current study, two metaheuristic approaches (i.e., GA and HS) were used to estimate the parameters of PDMs.

2.2.1. Genetic Algorithms

The GA is a stochastic search algorithm based on natural evolution and mechanisms of population genetics [8, 9]. The simple ideas of GA are based on the biological processes of survival and adaption. Only the best of the population is allowed to survive and propagate to successive generations. Therefore, the GA does not require derivatives of the objective function to solve complex and discontinuous optimization problems. A number of GAs are introduced, but the following general description encompasses most of the important features.

The analogy with nature is established by the creation of a set of solutions, called a population. The initial population of solutions is usually chosen at random and allowed to evolve over a number of generations. Each individual in a population is represented by a set of parameter values that completely describe a solution and undergo constant change by means of genetic operations of reproduction, crossover, and mutation. At each generation, the fitness of individuals with respect to the objective function is calculated for reproduction and propagated to the next generation. Based on the fitness, individuals with relatively high fitness (called parent) are selected for reproduction of the next generation. For example, as shown in Figure 1(a), there are two parents selected with binary cording. The strings, including the last three digits of two parents (i.e., 110 and 011), are recombined to produce offspring that will comprise the next generation. Then, the parents are replaced in the population by the offspring to keep a stable population size. The recombination operation is usually called the crossover. The offsprings (new generations) have a higher average fitness than their parents (previous generations). Occasionally, mutation is introduced into the population to prevent the convergence to a local optimum and help generate unexpected directions in the parameter space, as shown in Figure 1(b). A fixed fitness of generations has been created when a generation has reached the highest fitness and there is no further improvement with repeated iteration [23].

(a) Crossover

(b) Mutation

Figure 1 

Illustration of crossover and mutation.

The GA method has been widely applied in a variety of engineering optimization problems. It has been established that the GA approach is an attractive alternative to solve optimization problems with nondifferentiable, nonlinear, and multimodal objective functions [9, 25, 26].

2.2.2. Harmony Search

The HS developed by Geem et al. [24] is a phenomenon-mimicking algorithm inspired by the improvisational processes of musicians. The algorithm is based on natural musical performance processes in which a musician searches for a better state of harmony. Assume that the optimization problem is Minimize𝑂_𝑏(𝑎),Subjectto𝑎_𝑖∈Low_𝑖,Up_𝑖,𝑖=1,…,𝑁, (1) where 𝑂_𝑏(𝑎) is an objective function; 𝑎 is the set of each decision variable 𝑎_𝑖; 𝑖=1,…,𝑁 (representing probability parameter); 𝑁 is the number of decision variables (the number of parameters for a probability distribution); and Low_𝑖 and Up_𝑖 are the upper and lower limits, respectively, of the decision variable 𝑎_𝑖.

To solve the optimization problem, the procedure of the HS is as follows.(1)Generate the harmony memory (HM) randomly up to the harmony memory size (HMS) from a uniform distribution, denoted as 𝑎^𝑗_𝑖∼𝑈Low_𝑖,Up_𝑖, (2) where 𝑖=1,…,𝑁; 𝑗=1,…,HMS; and 𝑈[𝑎,𝑏] represents the uniform distribution ranging from 𝑎 to 𝑏. The generated HM is 𝑎HM=¹¹𝑎^1l₂⋯𝑎¹_𝑁𝑎²₁𝑎²₂⋯𝑎²_𝑁⋮𝑎^HMS₁𝑎^HMS₂⋯𝑎^HMS_𝑁=𝑎¹𝑎²⋮𝑎^HMS. (3)(2)Improvise a new set (𝑎_𝑖) by 𝑎_𝑖=𝑎_𝑖∈𝑎¹_𝑖,𝑎²_𝑖,…,𝑎^HMS_𝑖𝑎w.p.HMCR_𝑖∼𝑈Low_𝑖,Up_𝑖otherwise, (4) where HMCR is the harmony memory considering rate for all the variables 𝑎_𝑖 and 𝑖=1,…,𝑁. Equation (4) implies that each decision variable of a new harmony set is sampled from the same variable of HM in (3) for the probability of HMCR. Otherwise, generate the harmony set from the uniform distribution as in (2).(3)Adjust each variable of the improvised new set (𝑎_𝑖) with the probability of the pitch adjusting rate (PAR) as 𝑎^∗_𝑖=𝑎_𝑖𝑎+𝜀w.p.PAR_𝑖otherwise, (5) where 𝜀 is generated from 𝑈[−ℎ,ℎ]. Here, ℎ  is the arbitrary distance bandwidth. If ℎ is large, the variability of the adjusted value 𝑎^∗_𝑖 is also large.(4)Update the HM by replacing the worst harmony corresponding to the worst objective function with the improvised and adjusted new set.(5)Repeat steps (2) to (4) until the termination criterion has been met.

Following the empirically based HS parameter range is recommended to produce a sufficient solution of 0.7–0.95 for HMCR, 0.2–0.5 for PAR, and 10–50 for HMS [24].

2.2.3. Optimization Function

The derivative-based method usually uses the integral of the square of the error (ISE) as an optimization function to find a global solution because a derivative of ISE can be obtained relatively well. The metaheuristic approaches use a stochastic random search instead of a derivative search to find global optimization, so that various forms of the optimization function based on the integral of the absolute magnitude of the error (IAE), the integral of the time-absolute error (ITAE), and ISE are applied.

In this study, we use two metaheuristic approaches to estimate the parameters of PDMs. To attain the optimum result, an equivalent objective function (or target function) is used with two metaheuristic approaches. The objective function in this study is constructed using ISE with observed and estimated values and is expressed by Minimize𝑂_𝑏=^𝑛_𝑖=1𝑥_𝑖−𝑋_𝑖𝑥_𝑖², (6) where 𝑂_𝑏 is an objective function; 𝑥_𝑖 is the 𝑖th ordered observation value; and 𝑋_𝑖 is the estimated value from the corresponding 𝑖th empirical cumulative probability of the selected PDM.

2.3. Statistical Test Criteria

Three test criteria are used to assess the adequacy of the proposed parameter estimation method: the correlation coefficient (CC), coefficient of efficiency (CE), and root mean square error (RMSE) [27, 28]. The mathematical forms of these criteria are as follows: ∑CC=^𝑛𝑖=1𝑥_𝑖−𝑥𝑋_𝑖−𝑋∑^𝑛_𝑖=1𝑥_𝑖−𝑥²∑^𝑛_𝑖=1𝑋_𝑖−𝑋^21/2, (7) where 𝑥 and 𝑋 are the mean of the observed data and the mean computed data from the fitted PDM, respectively. The range of CC is −1≤CC≤1. Here, the CC is also equivalent to the Filliben Q-Q correlation test, which was proposed by Filliben [29] as a test of normality. Vogel [30] extended this approach to lognormal and Gumbel distributions for the goodness-of-fit test and GEV [31]: ∑CE=1−^𝑛𝑖=1𝑥_𝑖−𝑋_𝑖²∑^𝑛_𝑖=1𝑥_𝑖−𝑋²1,−∞≤CE≤1,RMSE=𝑛^𝑛_𝑖=1𝑥_𝑖−𝑋_𝑖^21/2,0≤RMSE≤∞. (8)

Note that higher CC and CE values and lower RMSE values represent better performance. Estimated statistics are described with boxplots. Boxes indicate the interquartile range (IQR), and whiskers extend to 1.5 IQR. The horizontal lines inside the boxes depict the median of the data.

3. Simulation Study

In the simulation study, we employ the GEV distribution model and estimate the parameters (i.e., scale (𝛼), shape (𝛽), and location (𝑥₀)) using the GA and HS techniques. For the simulation study, we produce 100 data sets using a GEV random number. An individual data set with a record length of 100 has the same parameters, which are 𝛼=30, 𝛽=0.1, and 𝑥₀=100. The strategies for the GA and HS are shown in Table 1.

Table 1 

Strategy for the GA and HS.

One hundred series are simulated for the parameter estimation of the GEV distribution model. The results of the parameter estimation based on the GA and HS methods are compared using boxplots, as shown in Figure 2. As shown in Figures 2(a) and 2(b), the 𝛼 and 𝛽 derived from the HS indicate a better performance than those of the GA. In addition, the variability of 𝛼 and 𝛽 derived from the GA is greater than that of the HS. However, as shown in Figure 2(c), the 𝑥₀ derived from the HS and GA shows similar results.

(a)

(b)

(c)

Figure 2 

Boxplots of estimated GEV parameters using 100 simulations based on GA and HS; (a) scale parameter, 𝛼; (b) shape parameter, 𝛽; and (c) location parameter, 𝑥0.

The statistical test criteria in (7) and (8) were computed using the estimated parameters. In Figure 3, the test criteria for the proposed metaheuristic approaches are shown. Recall that higher CC and CE values and lower RMSE values represent better performance. As shown in Figure 3(a), the CCs for the HS are higher, while the variability is narrower than the CCs for the GA. The CE shows a similar pattern to the CC, as shown in Figure 3(b). The RMSE for the HS shows lower values than that of the GA, as shown in Figure 3(c). The RMSE indicates that the HS results have more accurate values than the GA results.

(a)

(b)

(c)

Figure 3 

Boxplots of test criteria for estimated GEV parameters using 100 simulations based on GA and HS; (a) correlation coefficient; (b) coefficient of efficiency; and (c) root mean square error.

Although the metaheuristic approaches (e.g., the GA and ACO) are reasonable alternatives to solve hydrological problems, significant computation time is an obvious disadvantage [23]. Therefore, we investigate the computation time of two metaheuristic approaches. The computation time for GEV parameter estimation using 100 simulations based on the GA and HS is compared using boxplots, and the results are shown in Figure 4. The computation time for the GA is in the range of 10 to 20 minutes. However, the computation time for the HS is very short and is less than 1 minute. To determine computation time by varying conditions, the population size for the GA is changed to 100, 500, and 1000, while the harmony size for HS is also changed to 100, 500, and 1000. As shown in Figure 5, the computation time for the GA rapidly increases as the population size increases. Meanwhile, the fitness values of the objective function are only slightly improved. In contrast, the computation time for the HS rarely increases despite increasing harmony size, and the fitness values of the objective function remain consistent. Note that the fitness values of the objective function for HS with a small harmony size are better than the fitness values of the GA. The HS approach for estimation of the GEV parameter produces more reliable results than the GA approach and uses less computation time.

Figure 4 

Boxplot of computation time for GEV parameter estimation using 100 simulations based on GA and HS techniques.

Figure 5 

Comparison of computation time for GEV parameter estimation based on GA and HS methods.

In addition, three PDMs (i.e., LN2, GAM2, and GUM) are used in a simulation study. The same procedure as the simulation study using GEV is employed for the parameter estimation of each PDM. The target values for the parameters of the applied PDMs and the estimated parameters by the GA and HS methods are reported in Table 2. In addition, the results of the three test criteria and the computation time for each PDM are summarized in Table 3. The results of additional simulation studies show that the two metaheuristic approaches are appropriate methods for the parameter estimation of PDMs. Furthermore, the difference between the GA and HS methods is that the HS method requires less computation time than the GA method while still providing reliable results.

Table 2 

Target values for the parameters of employed PDMs and estimated parameters by the GA and HS.

Table 3 

Comparison of test criteria for the three PDMs.

4. Case Study

In the current case study, we carried out precipitation frequency analysis for design storm with annual maximum hourly rainfall data recorded at 74 rainfall gauges in South Korea. The annual maximum hourly rainfall data were extracted from the Korea Meteorological Administration website, http://www.kma.go.kr/. The record lengths of extracted rainfall data range from 20 to 100 years (average record length is approximately 40 years), and the locations of the 74 rainfall gauges are presented in Figure 6.

Figure 6 

The locations of the employed 74 rain gauges in Republic of Korea.

Four PDMs (i.e., LN2, GAM2, GEV, and GUM) were applied to the annual maximum hourly rainfall data. Parameters corresponding to the four PDMs were estimated using the three parameter estimation approaches: ML, GA, and HS. To compare the three parameter estimation approaches, the quantiles for different return periods (𝑇=10, 25, 50, and 100 years) were estimated at the 74 rain gauges in South Korea. The variability of the quantiles for each PDM is shown in Figure 7. Note that quantiles imply the statistically estimated design storm at a given return period.

(a)   LN2

(b) GAM2

(c) GEV

(d) GUM

Figure 7 

Quantiles using four PDMs at 74 stations in Republic of Korea. (a) LN2, (b) GAM2, (c) GEV, and (d) GUM.

As shown in Figure 7, the results of frequency analysis for the employed rainfall data show that the quantiles estimated by the two metaheuristic approaches present similar distributions, except for the quantiles of 𝑇=50, 100 for LN2, as shown in Figure 7(a). However, the quantiles based on the ML method show a slightly different distribution when compared with the two metaheuristic approaches.

Table 4 summarizes the results of test criteria derived from the three parameter estimation approaches. The CCs for ML, GA, and HS represent similar values, and the CE represents a similar pattern to that of CC. However, RMSE for GA and HS shows lower values than the results of ML. The RMSE indicates that the two metaheuristic approaches for the parameter estimation of PDMs have a more accurate performance than ML.

Table 4 

Comparison of test criteria for the annual maximum rainfall data for the 74 rain gauges stations in Republic of Korea.

In addition, the results of parameter estimation for Seoul, Daejeon, Daegu, and Busan, which are four major cities of South Korea, are summarized in Table 5. Furthermore, values of the quantiles for the four cities are estimated, and the results are summarized in Table 6.

Table 5 

The results of the parameter estimation for the four major cities in Republic of Korea.

Table 6 

Quantiles for different return periods (𝑇=10, 25, 50, and 100 years) of the four major cities in Republic of Korea.

5. Summary and Conclusion

The HS was developed by Geem et al. [24] and has been applied to solve a variety of engineering problems in a number of previous studies. Because the HS adopts a stochastic random search to find the optimized best solution, initial value settings of decision variables and derivative information to reach global optimum are not required. Furthermore, after considering all of the existing vectors based on the HMCR and PAR, the HS generates a new vector, whereas the GA only considers two vectors. These features increase the flexibility of the HS and produce a better solution. Therefore, we applied the HS method to enhance design storm estimation in the current study.

The results of the simulation study based on the four PDMs (i.e., LN2, GAM2, GEV, and GUM) showed that the HS approach for the parameter estimation of PDMs produced reliable results when measuring test criteria (i.e., CC, CE, and RMSE). Furthermore, the fitness values of the objective function for HS were better than the fitness values of the GA with a small harmony size. Additionally, the computation time for the HS approach was less than that of the GA method.

Precipitation frequency analysis for design storm was conducted to assess the performance of the proposed method with the four PDMs (i.e., LN2, GAM2, GEV, and GUM) and the annual maximum hourly rainfall data recorded at 74 rainfall gauges in South Korea. The results of precipitation frequency analysis were compared with those of the ML and GA methods. The results in this study suggested that performances of the GA and HS were similar and presented more accurate results than the ML in precipitation frequency analysis for annual maximum hourly rainfall data. Accordingly, we conclude that the proposed parameter estimation method based on the HS approach is a useful alternative for the parameter estimation of PDMs, particularly when conventional methods cannot be applied to estimate the parameters of PDMs.

Appendices

Probability Distribution Models

A. Two-Parameter Lognormal (LN2) Distribution

The transformation of random variable 𝑋 is 𝑌=log_𝑎(𝑋), (A.1) where 𝑎 is the base of the logarithm.

Assuming that the mean and variance of 𝑌 are 𝜇_𝑦 and 𝜎²_𝑦, respectively, then 𝑓(𝑥)=𝑘√2𝜋𝜎_𝑦−1exp2log_𝑎(𝑥)−𝜇_𝑦𝜎_𝑦² (A.2) is the probability density function (PDF, i.e., 𝑓_𝑋(𝑥)=𝑑𝐹_𝑋(𝑥)/𝑑𝑥) of the LN2. The cumulative distribution function (CDF, i.e., 𝐹(𝑥)=𝑃(𝑋≤𝑥)) of the LN2 is 𝐹_𝑋(𝑥)=^𝑦_−∞𝑘√2𝜋𝜎_𝑦−1exp2log_𝑎(𝑥)−𝜇_𝑦𝜎_𝑦²𝑑𝑥. (A.3)

Note that 𝑘=1 if 𝑎=𝑒 (base-𝑒 logarithm) and 𝑘=log₁₀(𝑒) if 𝑎=10 (base-10 logarithm). In relation to the random variable 𝑋, 𝜇_𝑦 controls the scale and is called the scale parameter, while 𝜎_𝑦 controls the skewness and is regarded as a shape parameter. Additionally, note that in relation to the variable 𝑌, 𝜇_𝑦 is the location parameter, and 𝜎_𝑦 is the scale parameter.

The quantile function for LN2 corresponding to the nonexceedance probability 𝑞 is 𝑥_𝑇=exp_𝑎𝜇_𝑦+𝑧𝜎_𝑦, (A.4) where 𝑥_𝑇 is a quantile of return period 𝑇 for LN2 and 𝑧 is the standard normal variate corresponding to the non-exceedance probability 𝑞.

B. Two-Parameter Gamma (GAM2) Distribution

The PDF of the GAM2 is 𝑓(𝑥)=1|𝛼|Γ𝛽𝑥𝛼^𝛽−1−𝑥exp𝛼,Γ𝛽=^∞₀𝑧^𝛽−1𝑒^−𝑧𝑑𝑧, (B.1) where 0≤𝑥<∞ for 𝛼>0 and −∞<𝑥≤0 for 𝛼<0. The parameters 𝛼 and 𝛽 are the scale and shape parameters. The shape parameter 𝛽 is restricted to 𝛽>0, while 𝛼 may be positive or negative. Γ(𝛽) is the complete gamma function, which is the integral function. From (B.1), the CDF of GAM2 is derived to be 𝐹(𝑥)=1=𝑃𝛽,𝑧𝛽×𝛼Γ_𝑥0𝑥𝛼^𝛽−1−𝑥exp𝛼𝐹(𝑥)1𝑑𝑧,for𝛼>0,=𝑃𝛽,𝑧=1−𝛽×𝛼Γ_𝑥0𝑥𝛼^𝛽−1−𝑥exp𝛼𝑑𝑧,for𝛼<0. (B.2)

The function 𝑃(𝛽,𝑧) in (B.2) is the incomplete gamma function. Therefore, the quantile for GAM2 must be estimated numerically. The quantile 𝑥_𝑇 of return period 𝑇 for GAM2 corresponding to the nonexceedance probability 𝑞 is alternatively estimated using the following: 𝑥_𝑇=𝜇_𝑥+𝐾_𝑇𝜎_𝑥, (B.3) where 𝜇_𝑥 and 𝜎_𝑥 are the mean and standard deviation, respectively. 𝐾_𝑇 is the frequency factor, which is a function of return period 𝑇 and the skewness coefficient [2].

C. Generalized Extreme Value (GEV) and Gumbel (GUM) Distribution

The PDF and CDF of the GEV for a random variable 𝑥 are expressed as follows: 𝑓(𝑥)=1𝛼1−𝛽𝑥−𝑥₀𝛼^(1/𝛽)−1𝐹(𝑥)𝐹(𝑥)=exp−1−𝛽𝑥−𝑥₀𝛼^1/𝛽, (C.1) where 𝛼, 𝛽, and 𝑥₀ are the scale, shape, and location parameter, respectively. 𝛽 plays an important role such that if 𝛽=0, the distribution tends to resemble a type-1 or Gumbel distribution; if 𝛽<0, the resulting distribution is type-2 or Log-Gumbel distribution; and if 𝛽>0, it is a type-3 distribution or Weibull distribution. The quantile functions for GEV and Gumbel corresponding to the non-exceedance probability 𝑞 are given in the following: 𝑥_𝑇=𝑥₀+𝛼𝛽11−−𝑙𝑛1−𝑇^𝛽𝑥forGEV,_𝑇=𝑥₀1−𝛼ln−ln1−𝑇forGumbel, (C.2) where 𝑥_𝑇 is a quantile of return period 𝑇 for GEV and Gumbel.

Acknowledgment

The authors acknowledge that this work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean Government (MEST) (2013-0362).