Comparison of Performance between Genetic Algorithm and Sce-ua for Calibration of Scs-cn Surface Runoff Simulation

Global optimization methods linked with simulation models are widely used for automated calibration and serve as useful tools for searching for cost-effective alternatives for environmental management. A genetic algorithm (GA) and shuffled complex evolution (SCE-UA) algorithm were linked with the Long-Term Hydrologic Impact Assessment (L-THIA) model, which employs the curve number (SCS-CN) method. The performance of the two optimization methods was compared by automatically calibrating L-THIA for monthly runoff from 10 watersheds in Indiana. The selected watershed areas ranged from 32.7 to 5844.1 km 2. The SCS-CN values and total five-day rainfall for adjustment were optimized, and the objective function used was the Nash-Sutcliffe value (NS value). The GA method rapidly reached the optimal space until the 10th generating population (generation), and after the 10th generation solutions increased dispersion around the optimal space, called a cross hair pattern, because of mutation rate increase. The number of looping executions influenced the performance of model calibration for the SCE-UA and GA method. The GA method performed better for the case of fewer loop executions than the SCE-UA method. For most watersheds, calibration performance using GA was better than for SCE-UA until the 50th generation when the number of model loop executions was 3434 around 5150 (one generation has 100 individuals). However, after the 50th generation of the GA method, the SCE-UA method performed better for calibrating monthly runoff compared to the GA method. Optimized SCS-CN values for primary land use types were nearly the same for the two methods, but those for minor land use types and total five-day rainfall for AMC adjustment were somewhat different because those parameters did not significantly influence calculation of the objective function. The GA method is recommended for cases when model simulation takes a long time and the model user does not have sufficient time for an optimization program to search for the best values of calibration parameters. For other cases, the SCE-UA program is recommended for automatic model calibration.


Introduction
Recently, optimization techniques have become widely used in computer modeling applications.A major advantage of the application of optimization techniques with models is to obtain better results without significant effort and time.Optimization techniques can be divided into two groups: global search methods such as Genetic Algorithm (GA) [1] and Shuffled Complex Evolution (SCE-UA) [2] algorithm, and local search methods, such as the downhill simplex method [3], the pattern search method [4], and the rotating direction method [5].Numerous studies have reported that global search methods performed better than local search methods [6][7][8][9][10].
Numerous studies reported that the SCE-UA algorithm provided better estimates of the solution compared with other global and local search procedures for calibrating watershed models.Cooper et al. [6] evaluated three global optimization methods including the SCE-UA, GA and simulated annealing (SA).They found that the performance of optimization was influenced by the objective function and starting position of the optimization search, and the SCE-UA method provided better results compared with other methods.Franchini et al. [8] reported that the solutions from SCE-UA were stable but characterized by many parameter values set at the boundary of their potential range.Kuczera [11] indicated the reason the SCE-UA is better than GA due to the periodic global sharing of information between the local simplex searches.Thyer et al. [10] revealed that SCE-UA's use of multiple complexes and shuffling provided a more effective search of the parameter space than GA's single simplex.
Although most optimization applications for automated model calibration continue to be studied, recent research activity in optimization for environmental modeling is centered on systematic analysis and selection of the best alternative for environmental management strategies.Vasques et al. [12] applied GA with the first-order reliability method (FORM) for estimating the probability of system failure under a given wasteload allocation to achieve wasteload allocation solutions between treatment cost and reliability.Ngo et al. [13] proposed to optimize control strategies for reservoir operation by applying SCE-UA methods.Kerachian and Karamouz [14] identified operating rules for water quality management in reservoir-river systems using a stochastic GA based conflict resolution technique.
Park et al. [15] proposed application of an optimization method for effective water quality monitoring in a large river system.
New optimization algorithms have been developed and some algorithms demonstrated better capability than GA or SCE-UA for some problems [16,17].However, GA and SCE-UA algorithms are applied for automatic calibration and determination of alternative strategies for watershed management, while other researchers have enhanced algorithm capabilities.Table 1 summarizes case studies of GA and SCE-UA application over the last five years (2010-2014) .SCE-UA has been primarily applied for automatic calibration, while GA has been more widely used for decision-making associated with watershed management strategy.Therefore, GA and SCE-UA can be expected to be widely used for watershed and water quality management for selection of the best alternatives for management as well as automated model calibration.In this study, GA and SCE-UA methods coupled with the Long-Term Hydrologic Impact Assessment (L-THIA) model, which employs the curve number (CN) method, were applied to 10 watersheds.The two optimization methods were characterized by analyzing how they generated solution populations and found optimal values according to the number of iteration of model runs.This study might help for user to select between GA and SCE-UA or for developer to improve the algorithm.

Genetic Algorithm (GA)
The GA is based on the principle of natural selection and natural genetics combining an artificial survival of the fittest with genetic operators abstracted from nature [1].The GA method employs three distinct operations including reproduction, crossover, and mutation, which make the initial population to be evolved the next generation [64].The flow chart of the PIKAIA version 1.0, which is a self-contained, genetic-algorithm-based optimization subroutine developed at the High Altitude Observatory [65], is shown in Figure 1.The initial population of chromosomes is generated randomly, and the objective functions for each individual are calculated.Both parents in one iteration of the reproductive cycle are selected by PIKAIA's stochastic sampling mechanism for new individual generation.A chromosome-like structure for each selected parent is produced for use in crossover and mutation processes.One-point crossover between a pair of parent-individuals is implemented to create a pair of offspring-individuals.One-point crossover selects randomly point C in length L and exchanges the codes of the pairs as shown in Figure 2. Mutation is implemented here by exchanging the values of some randomly chosen grid cells within one parent.The most-fit individuals are selected to mate and reproduce.The procedures from generating an initial population to select to mate and reproduce are repeated until a specific condition is met.Please refer the user guide of PIKAIA 1.0 for more detailed information [65].In this procedure, the population of chromosomes or individuals is a group of generated values for an optimized parameter, and parents are selected values within the group based on the objective function in the hydrologic model.

Shuffled Complex Evolution (SCE-UA)
The SCE-UA was developed in the Department of Hydrology and Water Resources of the University of Arizona and combined the simplex procedure with the concept of controlled random search [3], competitive evolution [1], and a complex shuffling concept [2].The SCE-UA flow chart is illustrated in Figure 3. Number of complexes (p) and points (m) in each complex are initialized and sample size is s = p × m. s points are generated in parameter spaces randomly; function values are calculated at each s points; and s points are ranked from the worst criteria value to the best criteria.The s points are partitioned into p complexes containing m points in partitioning into complexes.Each complex is evaluated according to the competitive complex evolution (CCE) algorithm outlined separately.The sample population is sorted, and the sample population is shuffled into p complexes.Convergence and the reduction in the number of complexes are checked.In the CCE algorithm, each complex is evaluated and is illustrated in

Long-Term Hydrologic Impact Assessment Model (L-THIA)
The L-THIA model was developed to evaluate the impact of land use change and long-term characteristics of watershed hydrology and nonpoint source pollution by combination of CN and event mean concentration (EMC) methods [66].The L-THIA model is available as a GIS version (Figure 5) and Web-L-THIA version available [67].The L-THIA model is primarily used in watershed management to simulate various hydrologic conditions for land use change [68][69][70][71] and loads of various nonpoint source pollution [72][73][74]; evaluate various scenarios for land use management [75,76]; develop total maximum daily loads (TMDLs) [77]; and evaluate low impact development (LID) [78][79][80].

Materials and Methods
SCE-UA and GA were linked to the L-THIA model for automated calibration of surface runoff.Ten watersheds in Indiana were selected and calibrated for surface runoff for comparison of the two methods.CN values and amount of total 5-day rainfall for antecedent moisture condition (AMC) adjustment were calibrated automatically using an approximately 10-year calibration period for each watershed.During the automated calibration process, the generated individuals and their objective values were extracted and compared to review how each optimization method searches for the optimal value by using various graphical techniques.Calibrated CN values and 5-day rainfall for AMC adjustment from the two methods were graphically plotted.

Modeling Approach
CN value is assigned to the combination of land use type and hydrologic soil group.For more exact simulation, spatial rainfall distribution needs to be considered by L-THA for runoff simulation.Thiessen polygons for rainfall gages, land use classification, and hydrologic soil group layers were overlapped and CN values were assigned for each Thiessen polygon area.The surface runoff for different Thiessen polygon areas were simulated separately and summed.The L-THIA model predicted runoff based on daily simulation by considering AMC classification and the objective function was calculated monthly because the L-THIA model does not consider routing processes since the main application of the L-THIA model is for long-term simulation.The Nash-Sutcliffe (NS) coefficient is widely used for evaluating hydrological simulation performance [81] and was used as the objective function in this study: where Qobs is the observed direct runoff; Qsim is the simulated direct runoff; and _ Q obs is average observed direct runoff.

Optimization Approach
FORTRAN codes of GA for PIKAIA version 1.0 and SCE-UA for version 2.2 were used to develop automated calibration of the L-THIA model.In each subroutine, which calculates the objective function of the two optimization programs, the surface runoff was simulated by an L-THIA algorithm, and the objective function value was calculated using simulated and observed runoff data.
AMC is divided as AMC I for dry conditions, AMC II for normal conditions, and AMC III for wet conditions (Table 2).CN values for AMC II and 5-day rainfall for AMC adjustment were calibrated.For calibrating CN values, CN values were changed by adopting multiplication factors for each land use type for default CN values to maintain the relationship between CN value for a given land use and its hydrologic soil group.Therefore, the multiplication factors were optimized, and calibrated CN values were calculated by multiplying the multiplication factors and default CN values (6).The default CN values and 5-day rainfall for AMC adjustment are shown in Tables 2 and 3.The range of optimized CN values was set as ±20% of default CN values to avoid searching extreme CN values.

= ×
(2 where CNcal is calibrated CN value; Px is a multiplication factor for each land use; and CNdef is default CN value.The range of total 5-day antecedent rainfall for adjusting CN for AMC was evaluated, and the range of rainfall amounts for the growing season was set higher than for the dormant season.CN value for each AMC is modified as follows: where, CNI is CN value for AMC I and CNIII for AMCIII. Once the optimization algorithm set the CN II value using Equation ( 2) and amount of total 5-day antecedent rainfall for AMC, L-THIA determines AMC by amount of total 5-day antecedent rainfall before a rain day, adjusts CN if AMC is I or III using Equation (3) or (4), and simulates surface runoff.
Generated individuals in both GA and SCE-UA were set as 20,000 to fairly compare performance of the two methods.Therefore, the number of individuals in a population and generation was set as 100 and 200 times for GA optimization, respectively, and number of complexes, points per complex, and maximum allowed trials were set at 23, 23, and 20,000 for SCE-UA optimization, respectively.

Study Area and Data Preparation
The L-THIA model linked with the two optimization methods was applied to 10 Indiana watersheds to compare optimization performance of the two methods (Figure 6).Brief descriptions for the 10 watersheds are provided in Table 4.The selected watershed areas ranged from 33 to 5811 km 2 .Five watersheds (WD#1 to WD#5) had crops as the primary land use, and forested area was dominant for watersheds WD#7 to WD#9.WD#6 and WD#10 had pasture/grass and urban as the primary land use, respectively.Major hydrologic soil groups for the study watersheds were B and C. The State Soil Geographic (STATSGO) database soil layer was obtained from USEPA [83].Although the Soil Survey Geographic (SSURGO) database contains higher resolution soil information, the search process and behavior to identify best values when using SSURGO or STATSGO for the two optimization methods are the same.Hydrologic soil group layer was generated by combination of soil layer and hydrologic soil group information for each soil type.Hydrologic soil group was divided in four groups: A for high infiltration, B for moderate infiltration, C for low infiltration, and D for very low infiltration.
The National Land Cover Data (NLCD 2001) for use as the land use layer was obtained from the Multi-Resolution Land Characteristics Consortium web site [84].The land use type was reclassified into seven categories for L-THIA model application: water, industrial/commercial, row crop, low-density residential, high density residential, grass/pasture, and forest.
The 10 watersheds were selected such that streamflow gage stations operated by U.S. Geological Survey (USGS) were located at the outlets of these watersheds.Streamflow needs to be separated for obtaining direct runoff for the automated calibration process because the CN method estimates direct runoff.The processes of obtaining USGS streamflow data and separating direct runoff were performed using the web-based WHAT program [25] with the digital BFLOW filter method [85].The 169 rainfall stations from 562 stations within Indiana operated currently by the National Climatic Data Center (NCDC) were selected, and spatial distribution of rainfall stations as Thiessen polygon GIS data was generated using their location information.NCDC rainfall stations for each watershed are listed in Table 4. Daily precipitation data location information for rainfall stations were obtained from the NCDC website [86].Periods of 10 years for calibration were selected by considering available rainfall and streamflow data for each watershed.

Generated Individual between GA and SCE-UA
Figure 7 shows the generated individual that was a calibration parameter (PC) to adjust CN values for crop area, which is the dominant land use type (80% of total area) in WD#1.The WD#1 showed the highest NS values, so it was assumed to be a good example for comparing search results between the two methods.Total number of model loop executions was 20,000 for SCE-UA and 20,300 for GA (100 individual per evolution and 200 evolutions).In early stages of evolution, the GA method was very fast and efficiently went to optimal space but after that point, generated individuals were somewhat dispersed compared with the SCE-UA method as shown in Figure 7.The reason the generated individuals are focused on the optimal space until early stages of evolution and then dispersed is well explained by PIKAIA's user manual [70].The initial, random population is distributed more or less uniformly in parameter space, but by the 10th evolution, nearly all individuals are concentrated on the central peak.After the 10th evolution, the mutation rate begins to increase.As the mutation rate is further increased, optimized parameters will produce a large horizontal or vertical jump in solution space, named the crosshair pattern.The use of elitism is essential in PIKAIA for improving the best solution.The best solution keeps improving after the 10th evolution, despite the fact that the spatial distribution of the population as a whole exhibits a greater spread about the central peaks in the later phase of the evolutionary run. Figure 8 represents how the GA method generates individuals in each evolution stage.The x and y are individuals for optimized parameters of industrial land use (31% of total area) and high residential area (28% of total area) in WD#10, respectively.Until the 10th evolution, generated individuals are more focused on optimal space as shown in Figure 8c and the crosshair pattern was developed (Figure 8d).
Overall, generated individuals by the SCE-UA were more concentrated in an optimal space compared with those by GA (Figure 7). Figure 9 shows the occurrence frequency of generated individuals for the optimized parameter of cropland (80% of total area) in WD#1.The two optimization methods show the peak in the same range, but frequency at that range for the SCE-UA was about two times higher than that for GA.Figures 7 and 9 illustrate the SCE-UA method focuses more on optimal space for generating individuals.

Performance of GA and SCE-UA for Model Calibration
The generated individuals for SCE-UA fall around a high NS value more often than the generated individuals for GA as shown in Figure 10.The number of model loops executed was 20,000 for SCE-UA and 20,300 (100 individual per generation and 200 generations) for GA.The individuals in x and y were optimized CN multiplication factors for high-density residential and industrial area, respectively.Total number of model loop executions for automatic calibration with GA or the SCE-UA method for model calibration impacted the results.Until 5150 model loop executions, the GA method showed better performance in L-THIA calibration than the SCE-UA method, and after 5150 times, SCE-UA was better than GA (Table 5).
The performance of GA and SCE-UA linked to L-THIA shows the two global optimizations were able to optimize CN monthly runoff well with most watersheds achieving NS values above 0.5 for 10 year calibration periods (Table 5 and Figure 11).The NS values for calibration ranged from 0.531 to 0.819 for GA and from 0.531 to 0.820 for the SCE-UA method.The filled circles and unfilled circles in Figure 12a show the 1:1 scatter plot of the calibrated CN values for all land uses and for the dominant land use, respectively, found in each watershed.The unfilled circles were more concentrated around the 1:1 line.However, the calibrated CN values for minor land uses in each watershed were unlikely to be as sensitive to calibration as the dominant land uses because they do not significantly influence calculation of the objective functions.Therefore, some filled circles are not as well centered on the 1:1 line.Jeon et al. [45] proposed a global calibration method which obtains one CN parameter set for the best fit to several calibration watersheds rather than an arithmetic average of CN parameters for each calibration watershed to regionalize CN parameters, because unrealistic CN parameters from minor land uses and hydrologic soil groups and realistic CN parameters from major land uses and hydrologic soil groups could provide results that equally regionalize CN values.Figure 12b shows the 1:1 scatter plot for optimized total 5-day rainfall for AMC adjustment, and optimized values fall further from the 1:1 line than CN values.
The parameters identified by SCE-UA for CN values of minor land use type and 5-day rainfall for AMC adjustments were sometimes differently optimized than those for GA as shown in Figure 12.The optimized CN values by SCE-UA for industrial/commercial areas, high and low density residential areas, and total 5-day rainfall were less variable and more concentrated around the average value (Figure 13).

Discussion
The NS value for the SCE-UA performance was a little higher than that for GA performance when the number of model loop executions by SCE-UA was large enough to search the optimized parameter, and in this study 20,000 times was enough (Figure 9).Cooper et al. [6] evaluated global optimization methods including GA, SCE-UA, and simulated annealing (SA) using a conceptual rainfall runoff model for eleven storms on a Caribbean island, and they also reported that the SCE-UA method provided the best approximations to the optimal solution compared with other methods; for the NS value, the SCE-UA method frequently searches for an exact solution but the GA method searches only near the optimal solution.The SCE-UA method concurrently searches for the optimal value within parameter space using divided complexes as shown in Figure 3. Cooper et al. [6] and Kuczera [11] reported that SCE-UA's robust capability comes from the periodic global sharing of information among the local simplex searches.
However, the number of model loop executions by the optimization program influenced which global optimization program between GA and SCE-UA was better for automatic calibration of SCS-CN model runoff in this study.The two methods have different ways to search for optimal values.The GA method was faster and efficiently reached the optimal space by the 10th generation, and after the 10th generation, generated individuals were somewhat dispersed to improve the optimization (Figures 6 and 7).Wang et al. [87] found variations in the calibration parameter that were almost same pattern found in this effort and reported in Figure 8 for GA.Table 4 illustrated that GA converges very quickly, and this has been reported by many researchers [88][89][90][91].Although the SCE-UA method results were somewhat dispersed in initial stages of optimization until around 5000 model executions compared to the GA method, after that point the generated individuals were more focused on the optimal space than those by the GA method (Figures 7-9).The differences in searching optimal space between GA and SCE-UA methods influenced the optimization performance between the two methods based on the total number of model loop executions (Table 5 and Figure 11).
L-THIA takes a relatively short time to be automatically calibrated when linked with GA or SCE-UA.Total automatic calibration time for WD#4, which is the largest watershed (5844.1 km 2 ) in this study, took 83 min (4980 s) for 20,000 model runs for SCE-UA and 84 min (5040 s) for 20,300 model runs for GA on a Windows XP machine with Intel ® Core2 Quad CPU 2.40 GHz and 3.25 GB RAM.In most cases the L-THIA user can select the best optimization method, which is SCE-UA in this study, since the time for calibration is relatively short.
The characteristics of the two algorithms examined in this study might be extended to other models when the most sensitive parameters to the model results are selected as an optimized parameter.The execution of more complex and data intensive hydrological models need significant execution times, so calibrating the complex model becomes an important problem when applying the model on large watershed with long time calibration periods [92].The number of calibration loop executions with SWAT is very critical so use of an optimization method that is more efficient with fewer calibration loop executions may be required for automatic calibration.
Most researchers emphasized SCE-UA was better than GA when applied to conceptual rainfall-runoff models, which are relatively simple and require short times to run [6,8,10,11].Van Griensven and Bauwens [93] developed ESWAT, which is SWAT linked with SCE-UA and applied to the Dender River (1384 km 2 ).ESWAT was ended after 8000 runs.If it takes 10 min for a SWAT run, calibration by ESWAT would take 55 days for 8000 runs in this application.
Therefore, a modeler must sometimes select either the GA or SCE-UA method for model calibration considering the available time for calibration.If a model takes a long time to run and the modeler doesn't have enough time for an optimization program to identify the best values, the GA method is recommended.Otherwise, if the modeler has sufficient time for the optimization approaches to run fully, the SCE-UA method is recommended for more accurate model calibration.In this study, NS value was used as the objective function because it has been widely used for hydrologic model calibration, although NS value has a bias of large stream flow events [94].The performance of the two optimization methods with other criteria used in objective functions, such as relatively error or root mean square error, should be evaluated in future studies.

Conclusions
In this study, the performance of the GA and SCE-UA optimization methods was linked to the L-THIA model, which employs CN and EMC techniques to estimate watershed runoff and nonpoint source pollutant loads, and optimization performances were compared and evaluated through application to 10 watersheds with a 10-year calibration period.The CN values and total five-day rainfall for AMC adjustment were calibrated automatically.The generated individuals by the SCE-UA method were more focused on optimal values as they evolved.The generated individuals by the GA method were generated toward optimal space until the 10th evolution, but after the 10th evolution they were dispersed because mutation rates are increased and a crosshair pattern appeared.Increasing the mutation rate is essential for GA with PIKAIA version to improve the search for the solution.
The two global optimization methods demonstrated acceptable performance for CN runoff simulation, and the calibrated direct runoff values were well matched with observed monthly runoff data.Significant performance differences between the GA and SCE-UA methods for calibration were not found, but the NS values for the SCE-UA are slightly higher than those for the GA method.However, the number of model loop executions for optimization influences the optimization performance between the two methods because the approach to searching optimal space is different.The GA method is recommended when the number of model loop executions for optimization is limited, otherwise SCE-UA's search method, which globally shares the information from each local simplex search, provided slightly better performance.The results provide important information to help model users select an algorithm for optimization of model parameters and to help developers improve algorithms.

Author Contributions
Ji-Hong Jeon wrote the paper and applied L-THIA linked with GA and SCE-UA; Chan-Gi Park reviewed the references and provided some useful comments; and Bernard A. Engel designed and directed this research.

Figure 4 .
In this procedure, s points are a group of generated values for an optimized parameter, and function values are calculated from the objective function between observations and hydrologic model simulation values.
H: 0, L: 0, O: 3, C: 21, P: 6, F: 66, W:3 A: 0, B: 38,C: 57, I is industrial/commercial; H is high density residential; L is low density residential; O is developed open space; C is crop; P is pasture and grass; F is forest; and W is water; Bold and underline indicates the dominant land use type; ** A means high infiltration (low runoff); B means moderate infiltration(moderate runoff); C means low infiltration (moderate to high runoff); and D means very low infiltration (high runoff); Bold and underlined indicates dominant hydrologic soil group.

Figure 7 .
Figure 7.Comparison of generated individuals during evolutionary run between GA and SCE-UA.Number of model runs were 20,000 for SCE-UA and 20,300 (100 individual/generation and 200 generations) for GA.The individual is an optimized parameter of crop land in WD#1 (80% of total area).(a) GA generated individuals; (b) SCE-UA generated individual.

Figure 8 .
Figure 8. Evolution of the population during the GA evolutionary run and development of the crosshair pattern on part (D).The x and y are individuals for optimized parameters of industrial land use (31% of total area) and high residential area (28% of total area) in WD#10, respectively.(a) Initial population; (b) 5th generation; (c) 10th generation; (d) 15th generation; (e) 100th generation; (f) 200th generation.

Figure 9 .
Figure 9. Histograms for generated individuals by GA and SCE-UA.Number of model runs were 20,000 for SCE-UA and 20,300 (100 individual/generation and 200 generations) for GA.(a) Histogram for GA; (b) Histogram for SCE-UA.

Figure 10 .
Figure 10.NS values and individuals for two optimized parameters, which are industrial land use (31% of total area) and high residential area (28% of total area) in WD#10, respectively.Number of model runs were 20,000 for SCE-UA and 20,300 (100 individual/generation and 200 generations) for GA.(a) NS values for GA; (b) NS values for SCE-UA.

Figure 11 .Figure 12 . 1 :Figure 13 .
Figure 11.The performance of GA and SCE-UA for model calibration according to the total number of model runs.

Table 1 .
Case studies of GA and SCE-UA algorithm application for the last five years (2010-2014).

Table 3 .
Default CN values in this study.
Notes: * HD residential is high density residential; ** LD residential is low density residential.

Table 5 .
Automatic calibration performance of L-THIA based on monthly simulation.