Inverse Modeling of Soil Hydraulic Parameters Based on a Hybrid of Vector-Evaluated Genetic Algorithm and Particle Swarm Optimization

The accurate estimation of soil hydraulic parameters (θs, α, n, and Ks) of the van Genuchten–Mualem model has attracted considerable attention. In this study, we proposed a new two-step inversion method, which first estimated the hydraulic parameter θs using objective function by the final water content, and subsequently estimated the soil hydraulic parameters α, n, and Ks, using a vector-evaluated genetic algorithm and particle swarm optimization (VEGA-PSO) method based on objective functions by cumulative infiltration and infiltration rate. The parameters were inversely estimated for four types of soils (sand, loam, silt, and clay) under an in silico experiment simulating the tension disc infiltration at three initial water content levels. The results indicated that the method is excellent and robust. Because the objective function had multilocal minima in a tiny range near the true values, inverse estimation of the hydraulic parameters was difficult; however, the estimated soil water retention curves and hydraulic conductivity curves were nearly identical to the true curves. In addition, the proposed method was able to estimate the hydraulic parameters accurately despite substantial measurement errors in initial water content, final water content, and cumulative infiltration, proving that the method was feasible and practical for field application.


Introduction
Soil hydraulic properties are essential characteristics for describing field-scale water flow and solute transport processes [1][2][3].Therefore, the accuracy and rapid estimation of soil hydraulic parameters are vital [4].Several direct methods have been developed to estimate the hydraulic properties of unsaturated soils [5][6][7].Unfortunately, direct methods have various limitations such as being time consuming and tedious, as well as being difficult for use in the field because of spatial variations in soil hydraulic properties [8][9][10].
The development of computer technology has led to an increased interest in indirect methods.Zachmann et al. [11,12] were the first to study the estimation of hydraulic properties using transient flow experiments in the laboratory.Since then, studies have developed inverse methods for estimating soil hydraulic parameters using in situ or laboratory experiments and silico simulation, such as one-step outflow and multistep outflow experiments [13][14][15], and evaporation experiments [16][17][18]; studies have proposed estimation devices, such as the tension disc infiltrometer [19,20], the modified cone penetrometer [21,22], and the field multiple extraction device [23].
The nonuniqueness and instability of parameter solutions in inverse methods are the key problems considered in most relevant studies.Most inverse methods use multitarget optimization algorithms, which utilize various data to improve the stability of solutions.Toorman et al. [19] evaluated inverse methods related to the one-step outflow method from the perspective of various objective functions, including combinations of matric potential, water content, and outflow data; specifically, they discovered that the matric potential data could greatly improve the sensitivity of parameter estimation.Dam et al. [24] concluded that, based on the multistep outflow method, the use of cumulative output alone can easily lead to a nonunique solution; however, with the addition of water content, a reliable parameter solution can be obtained.Eching and Hopmans [25] and Eching et al. [26] determined that the calibration of the parameters of the soil water retention curve was greatly improved using cumulative transient outflow data with simultaneously measured soil matric potential data; furthermore, they discovered that the addition of automatically measured soil pressure head data provided unique parameters in their multistep outflow experiment.Vereecken et al. [27] analyzed the effect of four parameters on the estimation of hydraulic functions based on multistep outflow experiments; they obtained an improved estimation with the addition of moisture retention characteristic data sets obtained directly from an outflow experiment.Vrugt et al. [28] proposed the parameter identification method based on the localization of information, which can be used in multistep outflow experiments to improve the uniqueness and stability of inverse parameters.
Tension disc infiltration has the advantages of short measurement time, wide availability of measured data sets, and easy application.Given these advantages, an in silico experiment was conducted to simulate tension disc infiltration.In this study, we first analyzed the synergistic effects of θ s , α, n, and K s by the objective functions of cumulative infiltration (Q), infiltration rate (v), and final water content (θ final ) using the van Genuchten-Mualem [29] model.Then, based on a hybrid vector-evaluated genetic algorithm (VEGA) [30] and particle swarm optimization (PSO) method [31], we proposed a new inverse method of soil hydraulic parameters named the "two-step method" under in silico experiments of tension disc infiltration, which first searches the hydraulic parameter θ s by the objective function of θ final , and then searches the soil hydraulic parameters α, n, and K s using the hybrid VEGA and PSO method by the objective functions of Q and v. Subsequently, because of the existence of local multiminimum domains of the objective functions in a tiny range of the true values, we analyzed the uniqueness and multiminimum values of the inverse method, as well as the stability and robustness of the inverse method.Finally, we analyzed the feasibility and stability of the proposed method in the presence of measurement errors.

Unsaturated Flow Governing Equation
The governing flow equation for radially symmetric isothermal Darcy flow in a variably saturated isotropic rigid porous medium, assuming that the air phase plays an insignificant role in the liquid flow process, is provided in the following modified form of the Richards equation [32]: where θ is the volumetric water content (cm 3 cm −3 ); h is the matric potential induced by capillary action (cm); K(h) is the hydraulic conductivity (cm min −1 ); z is the depth from the soil surface (cm) measured with positive values in a downward direction; r is the radial coordinate (cm); and t is time (min).

Initial and Boundary Condition
In Equation (1), we did not consider root water uptake by plant roots and assumed that the porous medium was isotropic; in addition, we assumed that the initial water content (θ initial ) and the initial matric potential (h i ) were the same in the vertical direction.The following initial and boundary equations (Equations ( 2)-( 5)) were used when Equation (1) was solved numerically.
where h i is the initial matric potential (cm), h 0 is the time-variable supply matric potential (cm); and r 0 is the disc radius (cm).Equation ( 1), subject to the abovementioned initial and boundary conditions, was solved using a quasi-three-dimensional finite element model, SWMS-2D, developed by Šim ůnek et al. [33].The numerical solution was based on the mass-conservative iterative scheme proposed by van Genuchten et al. [34].

van Genuchten-Mualem Model
The soil water retention curve (SWRC), h(θ), and soil water conductivity curve (SWCC), K(θ), are described by Equations ( 6)-( 8) [29]: where θ r and θ s are the residual and saturated water content levels (cm 3 cm −3 ), respectively; K s is the saturated hydraulic conductivity (cm min −1 ); α is an empirical parameter that is inversely related to the air-entry pressure value (cm −1 ); n is an empirical parameter related to the pore-size distribution; l is an empirical shape parameter, normally equal to 0.5; m = 1 − 1/n; and S e is the effective saturation.

Formulation of the Inverse Problem
As evidenced in the previous analysis, h(θ) (Equations ( 6) and ( 8)) and K(θ) (Equations ( 7) and ( 8)) are highly dependent on the effective saturation, S e .In addition, from Equation ( 8), the parameters of θ r and θ s have collinear relationships, thus simultaneous inverse estimation of θ r and θ s is not possible.Moreover, the value of θ r is relatively small with tiny variation range and can be either obtained experimentally (e.g., measuring the water content on very dry soil) or defined as the water content at a large negative value of the matric potential (e.g., the permanent wilting point (h = −15,000 cm)) [29].Besides, θ r can also be obtained from the pedotransfer functions by mechanical composition, bulk density, etc. [10].Therefore, in the following inverse procedure, we set the value of θ r as constant.
The soil water content (θ), cumulative infiltration (Q), infiltration rate (v), and matric potential (h) during the process of the infiltration can be used for the estimation of soil hydraulic parameters.However, measuring the h during infiltration is difficult, and thus, we calibrated the soil hydraulic parameters based on θ, Q, and v.
The parameters α, n, and K s in the van Genuchten-Mualem model clearly have a complicated nonlinear relationship, and parameters α and K s vary greatly, which makes them even more difficult to use in an inverse model.Table 1 was obtained from RETC software [34].Figure 1 illustrates the relationship between α-K s (1a) and α-n (1b).Notably, as α increases, the corresponding K s sharply increase, and n increases linearly; the correlation of these parameters should be considered in future research.The ratio α/K s was used in the development of the inverse model.In Table 1 and other references, the scope of soil hydraulic parameters (θ r , θ s , α, n, K s ) and α/K s were determined in Table 2 (set l as 0.5), which covers the vast majority of soil conditions.The objective function was defined as the root mean squared error (RMSE) divided by the standard deviation (σ) of the observed Q, v, or θ, as shown in Equation ( 9): where Y * represents specific measurements, such as Q, v, and θ, at time t i ; Y denotes the corresponding model prediction data under parameter vector β; β is the vector of optimized parameters (α, n, and K s ); and m represents the number of measurement sets (of Q, v, and θ).
The objective function ψ can be used to construct the response surface, by varying parameters α, n, θ s , and K s , which can be formulated with one of Q, v, or θ in VEGA, or a combination of them in PSO.
Most relevant studies have used various combinations of objective functions with Q, v, and θ to estimate soil hydraulic parameters; such combinations can convert multiple objective optimization problems into a single-objective optimization problem.However, determining the reasonable weight for the combination objective function is difficult; furthermore, it is likely to increase the complexity of the response surface, as well as the difficulty of searching for a global optimization solution.We adopted a hybrid VEGA and PSO method to solve these inverse problems.

The Genetic Algorithm and Particle Swarm Optimization
Genetic Algorithm (GA) is a global stochastic search technique developed by Holland [35].It begins with initializing a population of candidate solutions randomly sampled from a feasible solution space.Each individual is an entity with a characteristic chromosome, which is the main carrier of genetic material.The binary coding process, which is the process that is most commonly used in GAs, must be implemented from the start.It can realize chromosome mapping from phenotype to genotype.In each generation, individuals are sorted and selected according to values calculated by a fitness function.Then, the system executes combination crossover and variation on all members of the population, using the genetic operators of natural genetics.The best chromosome of the last population can be decoded and regarded as an approximate solution to a practical problem.This process causes the population to evolve in a manner similar to natural evolution, as the most recent generation is more adapted to the environment than the previous generation was.After initialization, the population evolves through the aforementioned steps and approximates the solution more closely with each iteration.This iterative evolution of the parameter solutions ends only when an optimal goal of model evaluations has been reached.
Particle Swarm Optimization (PSO) is an evolutionary optimization algorithm developed by Kennedy and Eberhart [31], which is derived from a study on bird foraging behavior.Initial particles are produced in the solution space.Each particle (P i ) in PSO has three characteristics: its position (X i ) represents a potentially optimal solution, its velocity (V i ) represents the direction and speed of the particle's movement, and its fitness value (F i ) (calculated by the fitness function) represents the relative merit of the particle.The particle moves in the solution space and updates the its location by tracking the current optimal position (P best ) and the current global optimal position (G best ), that is the fitness of the position calculated by one individual and by all particles.In each iteration, each particle is updated; thus, the P best and G best collected from the particles in each generation are updated.Each particle's speed can be adjusted dynamically according to its motion and the motions of the other particles, so that the best solution can be achieved in the optimal solution space.
In PSO, each particle (P i ) can be updated using the following equations: where w < 1 is the inertial weight; c 1 and c 2 are learning factors; and r 1 and r 2 are random numbers uniformly distributed in (0,1).The basic GA with a single population is powerful that its solutions are applicable in most situations.The multi-populations algorithm, in which the individuals can be exchanged through the various sub-groups, is closer to the nature of human evolution and often achieves a more desirable result.GAs have global search capabilities and robust performance, but their optimization searches may have limited efficiency.PSO searches run faster and more accurately than GA searches in local domains, but PSO is susceptible to unstable optimized solutions compared with GA.

Multiobjective Optimization Method Using Hybrid VEGA and PSO Algorithms
In general, inverse problems actually have multiple response variables, such as cumulative infiltration (Q), soil water content (θ), and infiltration rate (v).As these variables may conflict with each other, minimizing a single objective function often results in unacceptable outcomes for other objectives.One classical approach is to combine the individual objective functions into a single composite function.A weight can be assigned to each normalized objective function to solve the multiobjective optimization problem.Consider a decision maker who wishes to optimize multiple objectives, such that the objectives are incommensurable; the decision maker has no clear preference for one or the other.Therefore, in many practical problems, precisely and accurately selecting these weights can be difficult, even for someone familiar with the problem domain.
To solve the multiobjective problems, a VEGA [30] can provide a straightforward and efficient approach compared with the single-objective GA.In this algorithm, the whole population is divided into several subpopulations, which are evaluated with respect to the various objective functions.Normally when a VEGA operates, the whole population is randomly initialized and uniformly divided into K subpopulations of equal size.Each solution in a subpopulation is assigned one fitness value based on each objective function.Solutions are selected from these subpopulations using proportional selection.Crossover and mutation are performed on the new population in the same manner as for a single-objective GA.In addition, the lost individuals in the selection process are complemented by newly created random individuals and the elite population that remains from the previous generation.
Multiobjective functions can be allocated separately by VEGA; the optimization results always fall into the fittest Pareto sets but not the unique, precise solutions [36].In fact, VEGA can manage the sufficiency and equilibrium of all objective functions, but the shortcoming of this algorithm is that the solution tends to converge to extremes for each objective.Because of this, PSO is introduced to modify the extreme population towards more suitable individuals, but it cannot ensure whether the best solution is the global minimum.Therefore, the algorithms of VEGA and PSO should be rapidly operated alternately to maintain the balance between global searching and the local accuracy of minima.The detailed steps of the hybrid VEGA-PSO algorithm in this inverse problem are as follows: Step 1. Generate the initial population.The initial population is generated using a binary coding system with total number of K × N; K is the number of subpopulations classified by VEGA, while N is the number of individuals in each subpopulation.Step  One cycle of this hybrid VEGA-PSO algorithm (Figure 2) finishes when the iteration numbers m 1 in VEGA and m 2 in PSO have been reached.If the process reaches the defined number m 3 (Gen), the final termination condition is met and the final solution is output.

Data Generation
This study conducted an inverse estimation of the soil hydraulic parameters for four soils (sand, loam, silt, and clay from Table 1) under three initial levels of water content, θ initial , (S e = 20%, 40% and 60%).The infiltration data used in this study were generated using SWMS-2D software [33].Only the first three hours of the infiltration process were considered with the time-variable supply matric potential, h 0 , which are described in Equation (12).The initial h i were calculated by the corresponding parameters in Table 1. Figure 3 illustrates the results of Q under θ initial = 40% effective saturation degree.

Similarity Evaluation Criteria between Estimated and True SWRC and SWCC
To evaluate the similarity of the estimated (by inverse solution) and true soil water retention curves (SWRCs) and soil water conductivity curves (SWCCs), a certain amount of points were selected on the estimated (by inverse solution) and the true curves, with the same suction under SWRC and the same conductivity under SWCC.Subsequently, the soil water content levels of various points were calculated.Therefore, the similarity between the estimated and true SWRCs and SWCC were quantified using the following statistical indices: the root mean square error (RMSE), percent bias (PBIAS), and Nash-Sutcliffe coefficient (NS) [37].These indices are defined as follows: where m is the number of selected points in each SWRC and SWCC, which was set as 50 in our study; θ i est is the ith point of soil water content on the estimated SWRC and SWCC; and θ i true is the ith point of soil water content on the true SWRC and SWCC, with the same suction and conductivity.

Response Surface and Inverse Solutions
Based on the aforementioned analysis, we concentrated on the soil hydraulic parameters of θ s , α, n, and K s , with the θ r set as constant.To more clearly analyze the mutual interactions in objective function ψ of the four parameters θ s , α, n, and K s , the response surfaces for various parameter planes (α-n, α-K s , n-K s , n-θ s , α-θ s , and K s -θ s ) for four types of soil (sand, loam, silt, and clay) with three initial soil water content levels (S e = 20%, 40%, and 60%) were calculated.Each parameter domain was evenly divided into 50 discrete points, resulting in 2500 (50 × 50) grid points for each response surface.All operations were conducted using an Intel ® XEON ® CPU E5-2683 2.00 GHz processor and 32 GB of RAM in a Windows 7 Ultimate environment; a 28-core server was available for us to compute the algorithm.The objective function ψ(θ final ) was defined using the water content at the end of infiltration (θ final ). Figure 4 illustrates the response surfaces of parameter planes, including α-n (Figure 4a), α-K s (Figure 4b), n-K s (Figure 4c), n-θ s (Figure 4d), α-θ s (Figure 4e), and K s -θ s (Figure 4f), with the remaining two parameters set at the true values in the objective function of ψ(θ final ).For visualization, the logarithmic coordinates were selected for α and K s .From observing the first three contours in Figure 4a-c, the three planes α-n, α-K s , n-K s and the minimum values of ψ(θ final ) all clearly fall in a significant narrow and long valley.This means that even if the other two parameters are determined, the stable optimal values of α, n, and K s are still difficult to determine from ψ(θ final ).In Figure 4d-f, all three response surfaces, θ s -n, θ s -α, and θ s -K s , display smooth uniform rings and have obvious extreme points in each minimum circle.This feature suggests that it was possible to use the objective function, ψ(θ final ), to obtain the true value of θ s .

Analysis of the Objective Functions ψ(v) and ψ(Q)
If the value of θ s can be determined with objective function ψ(θ final ), similar to the abovementioned analysis, the mutual interactions of α-n, α-K s , and n-K s can be determined with objective functions ψ(v) and ψ(Q); the results of those calculations are illustrated in Figures 5 and 6.In Figure 5a-c, α-n, α-K s , and n-K s are depicted; all response surfaces evidently have global optima; therefore, after the value of θ s is determined, finding the true values of α, n, and K s is possible using the optimization method under the objective function ψ(v).Furthermore, the response surface ψ(v), with the parameters α, n, and K s , was re-estimated near its true value and the results are shown in Figure 5d-f.The response surfaces of α-n (Figure 5d) and n-K s (Figure 5f) can be observed to have obvious global optimal solutions, whereas for α-K s (Figure 5e), global optimal solutions exist in the response surfaces with numerous optimal bubbles in the valley.This creates difficulty in precisely searching for the true values of α, n, and K s only using ψ(v).
In Figure 6a-c, the response surfaces of α-n (Figure 6a), α-K s (Figure 6b), and n-K s (Figure 6c) clearly reveal that the minimum values of ψ(Q) were concentrated in a long narrow valley.Therefore, only using the objective function of ψ(Q) to search for the true values of α, n, and K s directly is difficult, but it can compress α-n, α-K s , and n-K s in a narrow space.Furthermore, similar to the aforementioned analysis, the response surface of ψ(Q) under the parameters α, n, and K s was recalculated to approximate its true value and the results are shown in Figure 6d-f.This further confirmed that the role of the objective function ψ(Q) is compressing α-n, α-K s , and n-K s in a narrow space in the search for α, n, and K s .Specifically, the response surface of the plane of α-K s (Figure 6e) indicated that the optimal value was decompressed in a distinct line-like interval, which was helpful for us in searching for the true values.
As a result, with synthetic analysis, the relationships between the parameter planes of α-n (Figure 5d) and n-K s (Figure 5f) with the objective function ψ(v), as well as α-K s (Figure 6e) with the objective function ψ(Q), determining the true values of α, n, and K s based on the multiple applications of the objective functions ψ(v) and ψ(Q) was possible.

The Procedure of Inverse Modeling
According to the aforementioned analysis, we first used the objective function of ψ(θ final ) to search for θ s with the GA.The values of θ final for the inverse modeling of typical soils with various θ initial values were simulated using SWMS-2D.The true parameters are shown in Table 1.For practical use, 24 grid points of θ final located at radial distance = 0, 5, 10, and 15, and depth = 2, 5, 10, 15, 20, and 35 were used for inverse parameters.According to the scope and features of these parameters, θ s and n were coded in arithmetic scaling and K s was coded in logarithmic scaling, whereas α was represented by the ratio value of α/K s , which was coded in arithmetic scaling.The population size (PopSize) of GA was 24 and the maximum number of iterations was 30.
Second, with the value of θ s set as in the abovementioned inverse model, the hybrid VEGA-PSO algorithm was used with two objective functions, ψ(v) and ψ(Q), comprehensively.The values of Q and v were selected every three minutes, with 60 points in total.The PopSize of PSO was 24, and each subgroup had two objective functions, ψ(v) and ψ(Q).The two subgroups of GA for ψ(v) and ψ(Q) had the same PopSize of 12.The coding mode of θ s , α, n, and K s were the same as the aforementioned step.For each individual, we set the generation iterations m 1 and m 2 as 3 in both GA and PSO for one coupling cycle, and the total number of coupling cycles was set at 30.The selection probability P s , the crossover probability P c , and the mutation probability P m were defined from GA optimization experience as the values 0.667, 0.80, and 0.20, respectively.In PSO, the inertial weight, w, and the learning factors, c 1 , and c 2 (in Equation ( 10)), were defined as the default values, equal to 0.729, 2.0, and 2.0 respectively.

Inverse Solution
To further verify the feasibility of the parameter inverse estimation approach, we conducted 12 cases including four types of typical soil from RETC (sand, loam, silt, and clay) under three initial water content levels (20%, 40%, and 60% degrees of effective saturation).Table 3 summarizes the results of the 12 cases of parameter inversion.Although small differences existed between the estimated and true values of each parameter, the estimates were close to the true values.

Analysis of the Inverse Solution
The true and the inversely estimated parameter values listed in Table 3 were used to draw the SWRC and SWCC, respectively, which are depicted in Figure 7. SWRCs are illustrated as a1, b1, c1, and d1; SWCC are illustrated as a2, b2, c2, and d2.Despite the small differences between the true and estimated values of α, n, and K s (Table 3), the true and the estimated SWRCs and SWCCs are almost identical, which indicates that the proposed inverse method is effective and robust.To evaluate the precision and robustness of the proposed inverse method, several performance indicators (RMSE, PBIAS, and NS) were adopted to estimate the differences between the true and estimated curves.Fifty points were selected from suction curves, h(θ), and conductivity curves, K(θ).Subsequently, the corresponding 50 soil water content levels were calculated.Table 4 illustrates the results of the error evaluation indicators, which indicate that our procedure is suitable for inverse modeling of soil hydraulic properties; this is because of extremely small MAE, RMSE, and PBIAS values for SWRCs (0.00039-0.000076cm 3 •cm −3 , 0.00071-0.00097cm 3 •cm −3 , and 0.0147-0.1603%,respectively) and also for conductivity curves (0.000004-0.00014cm 3 •cm −3 , 0.00003-0.00059cm 3 •cm −3 , and 0.2725-1.2131%,respectively), as well as a large NS (nearly 1.0).Furthermore, the inverse method is indicated to be able to estimate the SWRC and SWCC with precision and robustness.

Inverse Solutions with Measurement Errors
In real situations, some errors always exist, which are related to instrumentation, calibration, and other factors.Therefore, to assess the stability of the inverse solution, we considered a random error in cumulative infiltration and both the random error and system error for the initial and final water content levels, respectively.First, we set the random error as a 0.05 standard deviation (5.0%) and the two system errors as +0.02 and −0.02, respectively, for the θ initial and θ final .Subsequently, the two random errors were set as 0.02 and 0.05 of standard deviation (2.0% and 5.0%), and superimposed them to the cumulative infiltration.Similar to the error-free inverse method, all of the inversion computations were run again based on the proposed algorithm for the four representative soils (sand, loam, silt, and clay) with the θ initial (40% degree of effective saturation).

Inverse Solution Analysis Based on Initial and Final Water Content
Table 5 shows the inverse parameter values (θ s , α, n, K s ) of the four typical soils.The results suggested that superimposition of random and system errors on the error-free data result in only small deviations from the true parameters.Table 5 reveals that most values of the inverse parameters were even closer to their true values, which indicates that the proposed method was useful in practical application.Again, to assess the effective measurement error of the estimated values, the SWRCs and SWCCs were plotted based on the estimated and true values of the four typical soils in three error source cases, as shown in Figure 8 (initial water content) and Figure 9 (final water content).Both curves suggested that, although slight differences existed between the random and system error curves and the true curves, the effective measurement error was still acceptable.Table 6 illustrates several error evaluation indicators, including RMSE, PBIAS, and NS, to estimate the differences between inverse solutions with measurement errors and true values.Fifty points were averaged and selected from suction curves, h(θ), and conductivity curves, K(θ), and the corresponding 50 soil water content levels were calculated.Table 6 illustrates the results of error evaluation indicators, with extremely small RMSE (0.0018 to 0.03849 from h(θ), and 0.00039 to 0.02591 from K(θ)); small PBIAS (0.1401 to 14.1667 from h(θ), and 00847 to 10.0447 from K(θ)); and large NS (0.9308 to 0.9999 from h(θ), and 0.9319 to 1.0000 from K(θ)).The results indicated that the proposed inverse method is suitable for estimating soil hydraulic properties.

Inverse Solution Analysis Based on Cumulative Infiltration
In terms of cumulative infiltration, only random errors occurred; therefore, only two standard deviations of random errors, 0.02 and 0.05, were analyzed for the inverse solutions.Table 7 shows the results, which still suggest that small deviations existed between the estimated values with measurement errors and true values.Figure 10 depicts the SWRCs and SWCCs based on the inverse solutions within two standard deviations of random error.Small differences exist between the estimated values with measurement errors and the true values from these curves; their existence indicates that the inverse method was reasonable and robust for the inversion of soil hydraulic parameters.Furthermore, Table 8 demonstrates the values of RMSE, PBIAS, and NS between the true values and the estimated values with measurement errors.The small RMSE and PBIAS and large NS indicate the effectiveness and robustness of the inverse modeling.

Conclusions
In this study, we analyzed the inversion of soil hydraulic parameters for four typical soils (sand, loam, silt, and clay) under three initial water content levels (S e = 20%, 40%, and 60%).A new inverse method called the "two-step method" was proposed.According to the two-step method, the saturated water content, θ s , was primarily searched through GA with θ final .Subsequently, using the multiobjective optimization method based on the hybrid VEGA algorithm, the soil characteristic parameters α, n, and K s could be accurately estimated by cumulative infiltration and infiltration rate obtained by the SWMS-2D software for simulating water flow in two-dimensional variably unsaturated porous media.The results indicated that the proposed method is highly effective.In particular, compared with the traditional multitarget weighted sum optimization method, the hybrid VEGA-PSO algorithm method can overcome the difficulties of weight determination and search for the global optimal values efficiently and robustly.
Using the proposed method, accurate estimation of the hydraulic parameters remained difficult because the objective function had multilocal minimum values near the tiny range of true values.However, by comparing the SWRCs and SWCCs from the estimated and true values, we discovered that the curves were nearly identical.The results indicated that the proposed method can accurately estimate the SWRC and SWCC, although drawing a unique solution with the algorithm was difficult.Furthermore, considering that measurement errors of the initial water content, final water content, and cumulative infiltration inevitably exist in practical applications, comparison of the estimated SWRCs and SWCCs under measurement errors with the true curves proved that little difference existed.The results were acceptably accurate, which proved the proposed method to be robust and practical in the field.However, the inhomogeneity of soil textures and soil water content distribution require further investigation.

Figure 2 .
Figure 2. Flowchart of an entire optimizing cycle of the hybrid VEGA-PSO algorithm.

Table 1 .
Soil hydraulic parameters of the van Genuchten-Mualem model for 12 soils.
Note: Data for θ r , θ s , α, n, l, and K s were obtained from RETC.

Table 2 .
Inverse scope of soil hydraulic parameters.
2. Divide the population and calculate fitness in VEGA.VEGA divides the population into K subpopulations and the individuals in each subgroup are evaluated by corresponding fitness functions.Step 3. Sort and select separately.All the individuals are sorted and selected to eliminate unfit individuals based on the fitness values calculated in Step 2, but different fitness functions operate independently in each subpopulation.Step 4. Cross and mutate together.The populations selected in Step 3 are mixed together to execute the processes of crossing and mutation, in order to ensure the sufficiency and equilibrium of all fitness functions that VEGA includes.Step 5. Consider termination of VEGA.If the VEGA iteration number (Steps 3 and 4) is less than the specified value m 1 , then return to Step 3; otherwise, go to Step 6 (the PSO process).Step 6. Calculate fitness in PSO.All the individuals are evaluated in PSO by a composite fitness function with several measurement indices, to find minima.Step 7. Perform PSO operations.Each individual updates its velocity and position according to equations 10 and 11.Step 8. Consider termination of PSO.If the PSO process iteration number (Steps 6 and 7) is less than the specified value m 2, then return to Step 6; otherwise, return to Step 2 (VEGA process).

Table 3 .
Estimated parameters under three initial water content levels for four types of soil.

Table 4 .
Soil water content errors from estimated SWRCs and SWCCs for the four types of soil.

Table 5 .
Estimated parameters considering random and system errors on initial and final water content levels for the four typical soils.
Note: RE = random error; RE 5.0% means of standard deviation of random error equal to 0.05.SE = system error.

Table 6 .
Soil water content errors from estimated SWRCs and SWRCs considering both random and system errors on initial and final water content levels for the four typical soils.

Table 7 .
Estimated parameters considering random errors on cumulative infiltration for the four typical soils.

Table 8 .
Soil water content errors from estimated SWRCs and SWCCs considering random errors on cumulative infiltration for the four typical soils.: 2.0% and 5.0% means of standard deviation of random error equals to 0.02 and 0.05, respectively. Note