Particle Swarm Optimization and Shuffle Complex Evolution for Calibrating Xinanjiang Model Parameters

Xinanjiang model, a conceptual hydrological model, is well known and widely used in China since 1970s. Xinanjiang model consists of large number of parameters that cannot be directly obtained from measurable quantities of catchment characteristics, but only through model calibration. Parameter optimization is a significant but time-consuming process that is inherent in conceptual hydrological models representing rainfall–runoff processes. This study presents newly developed Particle Swarm Optimization (PSO) and compared with famous Shuffle Complex Evolution (SCE) to auto-calibrate Xinanjiang model parameters. The selected study area is Bedup Basin, located at Samarahan Division, Sarawak, Malaysia. Input data used for model calibration are daily rainfall data Year 2001, and validated with data Year 1990, 1992, 2000, 2002 and 2003. Simulation results are measured with Coefficient of Correlation (R) and Nash-Sutcliffe coefficient (E2). Results show that the performance of PSO is comparable with the famous SCE algorithm. For model calibration, the best R andE2 obtained are 0.775 and 0.715 respectively, compared to R=0.664 and E2=0.677 for SCE. For model validation, the average R=0.859 and average E2=0.892 are obtained for PSO, compared to average R=0.572 and average E2 =0.631 obtained for SCE.


INTRODUCTION
Xinanjiang model, a semi-distributed conceptual rainfall-runoff model was firstly developed in 1973 and published in English in 1980 (Zhao et al., 1980). Since its initial development in the 1970s, Xinanjiang model has been successfully and widely used in humid, semi-humid, even in dry area of China and elsewhere in the world for flood forecasting (Zhao, 1992;Zhao & Liu, 1995). According to Ren et al., (2006), Xinanjiang model has the following advantages: (a) The runoff depth are computed over partitioned sub-catchments: (b) evapotranspiration estimation by a three-layer method; (c) runoff components are separated into surface, subsurface, and groundwater flows, according to flow velocity; and (d) relations between partial model parameters on different time scales, can be conveniently transferred across temporal scales. One of the examples is deriving the outflow coefficients of the free-water storage to groundwater and subsurface flow in an hourly model, from a daily model based on the fact that daily hydrological data are more available than hourly data in the real world.
Xinanjiang model consists of large number of parameters that cannot be directly obtained from measurable quantities of catchment characteristics, but only through model calibration. According to Cheng et al, (2006) and Liu et al. (2009), calibration of the parameters is the main challenge in the development of hydrological models. In early days, the model calibration was performed manually, which is tedious and time consuming due to the subjectivities involved. Therefore, there is a need to employ automatic calibration techniques, which enables the hydrologist to rely less on subjective judgement for calibrating Xianjiang model.
In this study two Global Optimization methods named as Particle Swarm Optimization (PSO) and Shuffle Complex Evolution (SCE) were identified to calibrate Xinanjiang model parameters automatically. PSO was selected as Kuok (2010) found that Particle Swarm Method (PSO) is more reliable and promising to provide the best fit between the observed and simulated runoff for Hydrologic Tank model Calibration. Meanwhile, SCE method was selected as it was claimed that the most effective and efficient GOM for calibrating Hydrological Tank model parameters (Cooper et al., 1997;Chen et al., 2005).

PARTICLE SWARM OPTIMIZATION
PSO algorithm, a simple group-based stochastic optimization technique, was developed by Kennedy and Eberhart (1995). It is initialized with a group of random particles that were assigned with random positions and velocities. These particles learn over time in response to their own experience and the experience of the other particles in their group (Ferguson, 2004). Each particle keeps track of its best fitness position in hyperspace that has achieved so far (Eberhart and Shi, 2001) and accelerate towards its own personal best for every iteration. Then, the fitness value for each particle's is evaluated by calculating a new velocity term for each particle based on the distance from its personal best, and also its distance from the global best position.  "pbest" S e p t 25, 2 0 1 3 Once the particle achieved the best value, the particle stores the location of that value as -pbest‖ (particle best). The location of the best fitness value achieved by any particle during any iteration is stored as -gbest‖ (global best). The basic PSO procedure was shown in Fig. 1.
The particle velocity is calculated using Equation1.
Where Vi is current velocity,  is inertia weight, Vi-1 is previous velocity, presLocation is present location of the particle, prevLocation is previous location of the particle and rand() is a random number between (0, 1). c1 and c2 are acceleration constant for gbest and pbest respectively.  Once reached step 4 in SCE algorithm, Competitive Complex Evolution (CCE) will get into the procedure. CCE will reproduce new born and keep parents with stronger features, where each of the complexes will be developed separately before returning to the main stream of the procedure. All points in CCE have the chance to behave as parents in the way of reproduction. The duty of sub complexes is similar to parents except more than one parent can be utilized. Utilizing random method for the sub complexes enable users to conduct the search with complete converge of parameters domain. Schematic Diagram of SCE algorithm is presented in Fig. 2. Further details of SCE and CCE algorithm can refer to Duan et al. (1994).

XINANJIANG MODEL ALGORITHMS
The basic concept of runoff formation in Xinanjiang model is due to repletion of storage. Before the soil moisture content of the aeration zone reaching field capacity, runoff is not produced. Runoff is only generated at a point when the infiltration reached the soil moisture capacity (Zhao, 1983(Zhao, , 1992. Thereafter, the runoff is equal to excess rainfall without further loss. These processes are illustrated using tension water storage capacity curve, also named as parabolic curve of FC as presented in Equation 3 (Zhao et al., 1980).

(3)
Where WM'is the FC at a point that varies from zero to the maximum of the whole watershed WMM. LargerWM' means a larger soil moisture storage capacity in a local area and more difficult runoff generation. Parameter b represents the spatial heterogeneity of FC, as b =0 for uniform distribution and large b for significant spatial variation, and it is usually determined by model calibration.
According to Ren et al., (2006), three kinds of spatially heterogeneous distributions are taken into consideration in the Xinanjiang model and there are: (a) Uneven distribution of tension water storage capacity throughout the subcatchment is expressed by a parabolic curve for partial-area runoff generation; (b) Non-uniform distribution of free water storage capacity over partial area where runoff has been produced, is expressed also in terms of a parabolic curve for separation of runoff into surface flow, interflow and groundwater flow; c) The amount of free water storage is represented by the structure of linear reservoir for the sake of different velocities of different runoff components.
Meanwhile, the watershed-average soil moisture storage at time t ( , is the integral of between zero and , which is a critical FC at time t as presented in Equation 5 and Fig. 3: The critical FC ( ) corresponding to watershedaveragesoil moisture storage is presented in Equation 6.
When rainfall (Pt) exceeds evapotranspiration (Et), Pt is infiltrated into soil reservoir. Runoff (Rt) will only be produced when the soil reservoir is saturated (soil moisture reaches FC). As shown in Fig. 1, if the net rainfall amount (rainfall deduct actual evapotranspiration) in a time interval [t -1, t] is P-Et and initial watershed-average soil moisture(tension water) is Wt, then runoff yield in the time interval, Rt can be calculated as follows. S e p t 25, 2 0 1 3  Table 1. During the calibration, the parameter must satisfy the constraints of the Muskingum method for each channel of subbasin.  Exponential parameter with a single parabolic curve, which represents the non-uniformity of the spatial distribution of the soil moisture storage capacity over the catchment Im Percentage of impervious and saturated areas in the catchment Sm Areal mean free water capacity of the surface soil layer, which represents the maximum possible deficit of free water storage Ex Exponent of the free water capacity curve influencing the development of the saturated area Ki Outflow coefficients of the free water storage to interflow relationships Kg Outflow coefficients of the free water storage to groundwater relationships Ci Recession constants of the lower interflow storage Cg Recession constants of the groundwater storage

STUDY AREA
The selected study area is Bedup basin, with the approximate size of 47.5km 2 , located 80km from Kuching City, Sarawak, Malaysia. The catchment area is mainly covered with shrubs, low plant and forest, where there is no significant land use change over the past 30 years. Referring to topographic map produced by JUPEM (1975), the elevation of Bedup basin are varies from 8m to 686m above mean sea level. The length of Bedup river is approximately 10km. In general, Bedup basin is covered with clayey soils. Thus, most of the precipitation fails to infiltrate, runs over the soil surface and produce surface runoff. Part of Bedup basin is covered with coarse loamy soil, thus produced moderately low runoff potential.

Q=9.19( H ) 1.9 (7)
Where Q is the discharge (m 3 /s) and H is the stage discharge(m). These observed runoff data were used to compare the model runoff.
Daily areal rainfall data that obtained through Thiessen Polygon Analysis are fed into Xinanjiang model for model calibration and validation. Area weighted precipitation for BM, SN, SB, SM, ST are found to be 0.17, 0.16, 0.17, 0.18 and 0.32 respectively. Thereafter, the calibrated Xinanjiang model will carry out computation to simulate the daily discharge at Bedup outlet.

MODEL CALIBRATION AND VALIDATION
The model parameterization and model calibration is an iterative process. If the defined parameter selected produce poor calibration results, the model parameterization should be reconsidered by defining a simpler conceptual model that includes fewer calibration parameters. In contrast, if the model is not able to describe the response of the system sufficiently, key model parameters should be reconsidered and include other process descriptions in the calibration. Lower and upper limits on each parameter are usually used to define the parameter space. These limits are chosen according to physical and mathematical constraints, and information about physical characteristics of the system and from modeling experiences. The basic calibration procedure for Xinanjiang model using PSO and SCE algorithms is presented in Fig. 6. At the early stage of the calibration, the parameters of PSO and SCE that will affect the calibration results are pre-set. Various sets of daily rainfall-runoff data were calibrated to find the best model configuration. The objective function used is Root Mean Square Error (RMSE).
As the calibration process is going on, the initial parameters that set previously are changed to make the simulated runoff matching the observed one. The parameters investigated for PSO algorithm include:  Table 2. The accuracy of the simulation results are measured by Coefficient of Correlation (R) and Nash-sutcliffe coefficient (E 2 ). R and E 2 are measuring the overall differences between observed and simulated flow values. The closer R and E 2 to 1, the better the predictions are. The formulas of R and E 2 are presented in Equations8 and 9 respectively.
Where obs= observed value, pred = predicted value, = mean observed values and = mean predicted values.

PSO Algorithm Results
PSO algorithm achieved the optimal configuration at the RMSE of 2.3003. The optimal configuration for PSO algorithm was found to be 200 number of particles, max iteration of 150 and c1=1.8 and c2=1.8. The best R and E 2 obtained for calibration set are found to be 0.775 and 0.715 respectively as presented in Fig. 7. The 12 parameters of Xinanjiang model optimized by PSO algorithm can be found in Table 3.
The results show that runoff generated by Xinanjiang model optimized by PSO algorithm is controlled and dominant to 8 parameters include S, B, Im, Sm, Ex, Ki, Kg and Ci. In contrast, Dt, K, C and Cg are less sensitive to storm hydrograph generation. Fig.8 shows the validation results when the optimal configuration of Xinanjiang model optimized by PSO algorithm. As R is referred, the results obtained for Year 2000Year , 2003Year , 2002Year , 1992Year and 1990 are found to be 0.674, 0.649, 0.616, 0.616, 0.553 and 0.622 respectively. As E2 is used as level mark, the E2obtained are ranging from 0.550 to 0.623. The average R and E2 are yielding to 0.622 and 0.579 respectively.

SCE Algorithm Results
The optimal configuration of SCE algorithm is obtained at maxn=10000, kstop=100 and nsp1 of 101 with RMSE of 3.1536. The best R and E 2 are found to be 0.664 and 0.677 respectively for calibration set (as presented in Fig. 9). Optimal 12 parameters of Xinanjiang model obtained by SCE algorithm can be found in Table 4. The results indicate that Xinanjiang S e p t 25, 2 0 1 3 model optimized by SCE algorithm is dominant to 11 parameters except Dt. This implicate that all the 11 parameters include S, K, C, B, Im, Sm, Ex, Ki, Kg, Ci and Cg B and Cg are sensitive and have significant impact to runoff generation.