Estimating Soil Hydraulic Parameters during Ponding Inﬁltration Using a Hybrid Algorithm

: Accurate inversion of soil hydraulic parameters based on the van Genuchten–Mualem model has received much attention in soil science research. Herein, a hybrid algorithm method using particle swarm optimization and vector-evaluated genetic algorithm was used to invert the parameters θ s , α , n , and K s , with the objective functions of infiltration rate, cumulative infiltration, and soil water content. Then, numerical experiments were conducted on four typical soils at three initial water content levels (20, 40, and 60% effective saturation) to verify the accuracy of the inverse method. The results showed that the inversed soil water retention and conductivity curves were approximately the same as the real curves, with the root mean square errors of 0.00101–0.00192 cm 3 · cm − 3 , 0.00800–0.02519 cm 3 · cm − 3 , respectively, and both the Nash-Sutcliffe coefficients were approximately 1.0. Additionally, laboratory experiments were also performed to compare with the inversed parameters for verification, within small root mean squared errors and approximately 1.0 Nash–Sutcliffe coefficients. Furthermore, the method can also achieve acceptably accurate parameter inversion even with substantial measurement errors included in the cumulative infiltration, initial water content, and final water content. Thus, the method is effective and robust and found to be practical in field experiments.


Introduction
Numerical models describing the flow and transport of water and chemicals through the vadose zone have been increasingly and extensively applied in field-scale research. When these models are used to simulate saturated water flow and contaminant transport, it is critical to have a good understanding of unsaturated soil hydraulic properties, as the outcome from applying these models depends on the accuracy of model evaluation. Unsaturated hydraulic conductivity, K(θ) and soil water retention, h(θ) are often required to solve the Richards equation numerically. However, the extreme spatial heterogeneity of a subsurface environment confounds measurement of these hydraulic properties [1,2]. Hydraulic properties may exhibit significant variation over time due to changes in the ionic composition and concentration of a soil solution, the impact of soil crust and particle dispersion, shrink-swell phenomena in fine-textured soil, and cultivation or other agricultural activities [3].
The soil water diffusivity, D, and hydraulic conductivity, K, can be measured directly by laboratory and field methods as a function of water content or pressure head [4,5]. Direct measurement methods are relatively simple in concept, but several limitations restrict their practical usage [6]. As a result, the application of direct measurement methods in field gravity drainage experiments of layered profiles, or medium-and fine-textured soils is limited. In addition, methods that require equilibrium conditions, such as repeated steadystate flow situations, are time-consuming. Additionally, additional errors are introduced

Water Flow
Assuming that the air phase plays a non-significant role in the liquid flow process, the flow equation governing radially symmetric isothermal Darcy flow in a variably porous, rigid, isotropic, and saturated medium is given by the modified form of the Richards equation as follows: where t is time (min); r is the radial coordinate (cm); z is the depth of soil surface (cm) with a positive value indicating a downward direction; K(h) is the hydraulic conductivity (cm·min −1 ); h is the potential head (cm), and θ is the soil water content (cm 3 ·cm −3 ). In Equation (1), the root water uptake by plant roots was not considered and the porous medium was assumed to be isotropic. Additionally, it was also assumed that the metric head (hi) and initial water content (θi) were the same in the vertical direction.

Water Flow
Assuming that the air phase plays a non-significant role in the liquid flow process, the flow equation governing radially symmetric isothermal Darcy flow in a variably porous, rigid, isotropic, and saturated medium is given by the modified form of the Richards equation as follows: where t is time (min); r is the radial coordinate (cm); z is the depth of soil surface (cm) with a positive value indicating a downward direction; K(h) is the hydraulic conductivity (cm·min −1 ); h is the potential head (cm), and θ is the soil water content (cm 3 ·cm −3 ). In Equation (1), the root water uptake by plant roots was not considered and the porous medium was assumed to be isotropic. Additionally, it was also assumed that the metric head (h i ) and initial water content (θ i ) were the same in the vertical direction. Equation (1) was solved numerically, followed by the initial and boundary equations (Equations (2)-(5)): θ(r, z, t) = θ i , t = 0 (2) h(r, z, t) = h i , t = 0 h(r, z, t) = h 0 (t), 0 < r < r 0 , z = 0 (4) − ∂h(r, z, t) ∂z − 1 = 0, r > r 0 , z = 0 (5) h(r, z, t) = h i , r 2 + z 2 → ∞ where r 0 is the disc radius; h 0 is the time-variable supply pressure head, and h i is the initial pressure head. The SWMS [33] was used to solve Equation (1) under the previously mentioned initial and boundary conditions. Van Genuchten [31] proposed a mass-conservative iterative scheme on which the numerical solution was based.

Soil Hydraulic Properties
Prior to the numerical solution of the Richards equation, a parametric model of unsaturated soil hydraulic properties was selected by the parameter optimization approach. In this study, the unsaturated soil hydraulic functions of the van Genuchten-Mualem model [34] was adopted to describe the soil water conductivity curve (SWCC) and soil water retention curve (SWRC): where m = 1 − 1/n; K s is the saturated hydraulic conductivity (cm·min −1 ); l is an empirical shape parameter equal to 0.5; n is an empirical parameter associated with the pore-size distribution; α is an empirical parameter inversely associated with the air-entry pressure (cm −1 ); θ r is the residual water content level (cm 3 ·cm −3 ); θ s is the saturated water content level (cm 3 ·cm −3 ); S e is the effective saturation.

Inverse Problem and Hybrid Optimal Algorithm
As shown in Equations (7) and (8), h(θ) and K(θ) are highly related to the effective saturation, S e . Equation (7) demonstrates that the parameters θ r and θ s have a collinear relationship that cannot be inverted simultaneously. Moreover, θ r is relatively small and can be obtained through the application of a transfer function to other physical soil characteristics [35][36][37]. Hence, the value of θ r is set as the true value in the following inverse procedure. In the van Genuchten-Mualem model, the three parameters α, n, and K s exhibit a complex nonlinear relationship, especially the parameters α and K s which have a wide range of the order of magnitude, making it harder to be used in an inverse model.
During the infiltration process, potential head (h), infiltration rate (v), cumulative infiltration (Q), and soil water content (θ) can be used to estimate soil hydraulic parameters, but it is difficult to measure h. Thus, soil hydraulic parameters can be inverted based on θ, Q, and v in this study. Table 1 displays the soil hydraulic parameters of 12 typical soils obtained from RETC. Figure 2 details three kinds of parameter-pair relationships: α-K s (a), α-n (b), and v-K s (c). Figure 2a demonstrates that with the increase of α, the corresponding n increases linearly and K s increases sharply. Thus, it is necessary to consider correlations between these parameters in future research; otherwise, analysis of physical mechanisms is unreasonable and difficult to perform in the inverse model. Therefore, the ratio of α/K s Agronomy 2023, 13, 726 5 of 24 can be used in the inverse model. To improve the presentation of the inverse work, a comprehensive summary of the scopes of soil hydraulic parameters (K s , n, α, θ s , and θ r ) and the ratio of α/K s were mostly taken from the Unsated Soil Database (UNSODA) [38], as well as from a further 192 published available literatures, as shown in Table 2, covering the majority of soil conditions. A correct and relatively small range of parameters is critical for determining optimal values. Even when a restrictive range of α/Ks is employed, the range of Ks remains relatively large (6000 times of difference in magnitude), which is difficult to use in an optimal algorithm. During the infiltration process, the saturated conductivity (Ks) strongly influences the infiltration rate (v). The relationship between Ks and v for 12 typical soils with three depths (2, 3, and 5 cm) of ponding water is illustrated in Figure 2c. The infiltration rate was obtained using SWMS-2D, for which cumulative infiltration was simulated every 5 min. All the datasets were found to range from approximately 2 to 60, providing a narrow range (30 times), compared with the original range of Ks for the inverse algorithm. Hence, the ranges of both α/Ks and v/Ks were employed in the inverse procedure.    A correct and relatively small range of parameters is critical for determining optimal values. Even when a restrictive range of α/K s is employed, the range of K s remains relatively large (6000 times of difference in magnitude), which is difficult to use in an optimal algorithm. During the infiltration process, the saturated conductivity (K s ) strongly influences the infiltration rate (v). The relationship between K s and v for 12 typical soils with three depths (2, 3, and 5 cm) of ponding water is illustrated in Figure 2c. The infiltration rate was obtained using SWMS-2D, for which cumulative infiltration was simulated every 5 min. All the datasets were found to range from approximately 2 to 60, providing a narrow range (30 times), compared with the original range of K s for the inverse algorithm. Hence, the ranges of both α/K s and v/K s were employed in the inverse procedure.
During the parameter inversion process, the objective function ψ is minimized and used to structure a response surface. The function can be expressed by one or a combination of θ, v, and Q. Most relevant studies have formulated objective functions to estimate parameters using various combinations of θ, v, and Q, which can be used to transform into a single-objective optimization problem from a multi-objective optimization problem. However, it is difficult to determine a reasonable weight for the combination of objective functions, and a combined objective function may well enhance the difficulty of searching for a global optimization solution and the complexity of minimum response surface domains. A hybrid method of PSO and VEGA was adopted to solve the inverse estimation. The standard deviation (σ) of the observed data of θ, v, or Q dividing the root mean squared error (RMSE) was called the objective function: where m is the number of measurement sets (θ, v, and Q); β is the vector of optimized parameters (K s , n, and α); Y is the corresponding model prediction data under the parameter vector β; Y* is the specific type of measurement data at time t i . The hybrid optimal VEGA-PSO algorithm presents a highly effective optimization method. Holland [39] developed a global stochastic search technique called the genetic algorithm (GA). Eberhart and Kennedy [40] developed PSO, an evolutionary optimization algorithm from the research on bird foraging behavior. In this study, multiple objectives are still present in the inverse problems of real soil hydraulic parameters, including infiltration rate (v), soil water content (θ), and cumulative infiltration (Q). The main problem in multiobjective optimization is that inverting a single objective tends to cause unacceptable outcomes for other objectives as the objectives conflict with each other. It is often difficult to select weights precisely and accurately for the objectives in many practical problems, even for a researcher familiar with the problem domain. Compared to a single-objective GA, VEGA provides a straightforward and efficient way for solving multi-objective problems. More details about the hybrid optimal algorithm for inverse parameters may be found in our previous paper [32].
The above-mentioned analysis focused on the parameters of K s , n, α, and θ s , with θ r as its true value. The response surface of different parameter planes (K s -θ s , α-θ s , n-θ s , n-K s , α-K s , and α-n) was calculated for four types of soil (clay, silt, loam, and sand) with three initial soil water content levels (20, 40, and 60% effective saturation), to analyze the interaction in the objective function ψ of K s , n, α, and θ s . Each parameter domain was equally classified into 50 discrete points, and each response surface yielded 2500 (50 × 50) grid points. All operations were performed in a Windows 7 Ultimate environment using 32 GB of RAM and an Intel ® XEON ® CPU E5-2683 2.00 GHz processor. A 28-core server was used for computing through the algorithm.

Numerical Experiment
The SWMS-2D software [33] was used to generate infiltration data in this study. For the time-variable supply pressure head, h 0 (2 cm), the first 3 h (180 min) of the infiltration process were only taken into account. We concentrated on four soils: clay, silt, loam, and sand, with three θ initial , (represented by 20, 40, and 60% effective saturation) for the inverse modeling. The initial h i was calculated from each soil type's corresponding parameters, as shown in Table 1.

Laboratory Experiment
Experimental soil samples from depths of 0-60 cm were obtained from the Key Laboratory of Agricultural Soil and Water Engineering in Arid and Semiarid Areas of China (34 • 170 N and 108 • 040 E). Based on the USDA Soil Taxonomy System, soil samples were clay loam with particle size distributions of 27.30% 0-0.002 mm, 42.53% 0.002-0.020 mm, and 30.17% 0.02-2.00 mm. The experimental soils were air-dried, followed by screening through a mesh of 2 mm. Subsequently, they were compressed into a soil bin at a bulk density of 1.30 or 1.35 g cm −3 to simulate in situ bulk density. Water was mixed with the soil to achieve the desired experimental soil water conductivity (SWC) values of 10, 20, and 30% for the respective bulk density before the soil was compacted into the bin. A homogeneous soil profile was then created by loading and compacting the soil into the bin at a layer of 5 cm. The experiment was designed to determine the irrigation volume. In this way, a Marriotte bottle was placed to maintain an infiltration depth of 3 cm (Figure 3). Cumulative infiltration was recorded for each minute of infiltration. Eventually, soil samples were gathered from side holes, and SWC was determined by recording the weight loss of samples after 24 h of oven drying at 105 • C.

Laboratory Experiment
Experimental soil samples from depths of 0-60 cm were obtained from the Key Laboratory of Agricultural Soil and Water Engineering in Arid and Semiarid Areas of China (34°170 N and 108°040 E). Based on the USDA Soil Taxonomy System, soil samples were clay loam with particle size distributions of 27.30% 0-0.002 mm, 42.53% 0.002-0.020 mm, and 30.17% 0.02-2.00 mm. The experimental soils were air-dried, followed by screening through a mesh of 2 mm. Subsequently, they were compressed into a soil bin at a bulk density of 1.30 or 1.35 g cm −3 to simulate in situ bulk density. Water was mixed with the soil to achieve the desired experimental soil water conductivity (SWC) values of 10, 20, and 30% for the respective bulk density before the soil was compacted into the bin. A homogeneous soil profile was then created by loading and compacting the soil into the bin at a layer of 5 cm. The experiment was designed to determine the irrigation volume. In this way, a Marriotte bottle was placed to maintain an infiltration depth of 3 cm ( Figure  3). Cumulative infiltration was recorded for each minute of infiltration. Eventually, soil samples were gathered from side holes, and SWC was determined by recording the weight loss of samples after 24 h of oven drying at 105 °C.

Evaluation Criteria
To evaluate the similarity of the inverted and real SWRCs and SWCCs, a number of points (i.e., 50) were selected on the inverted and real curves. For SWRCs, points with the same suction were selected; for SWCCs, points with the same conductivity were selected. Afterwards, the soil water content levels were calculated at each point. The statistical indices of the Nash-Sutcliffe coefficient (NS), percent bias (PBIAS), RMSE, and mean absolute error (MAE) were used to quantify the similarity between the inverted and real SWRCs and SWCCs, which are defined as follows:

Evaluation Criteria
To evaluate the similarity of the inverted and real SWRCs and SWCCs, a number of points (i.e., 50) were selected on the inverted and real curves. For SWRCs, points with the same suction were selected; for SWCCs, points with the same conductivity were selected. Afterwards, the soil water content levels were calculated at each point. The statistical indices of the Nash-Sutcliffe coefficient (NS), percent bias (PBIAS), RMSE, and mean absolute error (MAE) were used to quantify the similarity between the inverted and real SWRCs and SWCCs, which are defined as follows: where θ i real is the ith point of the soil water content of the real SWRCs and SWCCs, with the same conductivity and suction; θ i inv is the ith point of the soil water content of the inverted SWRCs and SWCCs; m is the number (50) of selected points in the inverted and real SWRCs and SWCCs.

Response Surface and Uniqueness
The objective function surfaces of 12 typical soils had similar shapes (Table 1). Therefore, only the response surface of loam with 40% initial soil water content was shown in this paper. Figure 4 presents the response surface of parameter planes α-n (a), α-K s (b), n-K s (c), n-θ s (d), α-θ s (e), and K s -θ s (f), with the remaining two parameters as the real values, which were determined using the objective function ψ(θ final ). Clear visual observation was facilitated by selecting logarithmic coordinates for K s and α. Figure 4a shows a surface with a wide valley resulting from high n and low α, whereas Figure 4b,c show narrow and long valleys for α-K s and n-K s . As can be observed from the first three valley contours in Figure 4a-c, it remains difficult to determine the optimal stability values of K s , n, and α from ψ(θ final ) even though the other two parameters are identified. In Figure 4d-f, however, smooth uniform rings with distinct extreme values in each minimum circle were shown on all three response surfaces (θ s -K s , θ s -α, and θ s -n), suggesting that the true value of θ s may be obtained using the objective function ψ(θ final ).
The feasibility of applying ψ(Q) and ψ(v) as the objective functions to inverse the remaining parameters (K s , n, and α) through the response surface was then analyzed based on the present method. The results are shown in Figures 5 and 6. It can be seen from Figure 5a-c that the response surfaces corresponding to α-n, α-K s , and n-K s have obvious global optimal characteristics, indicating that it is possible to invert α, n, and K s according to the objective function of ψ(v) after θ s is determined. Furthermore, the regions where the optimized results (α, n and K s ) are close to the true values in each response surface (Figure 5a-c) are enlarged, as indicated in Figure 5d-f. It can be seen from Figure 5e, f that there are obvious global optimal solutions in response surfaces of n-K s and α-n. However, the response surfaces show some bubbles in the global optimal solution area, which indicates that the optimal solutions are not unique. Thus, it is demonstrated that the solutions are non-unique if only ψ(v) is applied as the objective function to inverse α, n, and K s simultaneously.
Moreover, taking ψ(Q) as the objective function, Figure 6a-c illustrates that the optimal solutions of n-K s , α-K s , and α-n response surfaces are narrow and long regions, which are relatively small compared to Figure 4a-c. Meanwhile, the above similar approach is adopted to enlarge the regions where the inversion results of α, n, and K s are close to the true values in each response surface (Figure 6a-c). It is further proved that the optimal solution range of the response surfaces of n-K s , α-K s , and α-n is smaller with ψ(Q) as the objective function, as shown in Figure 6d-f.
In summary, with the ψ(θ final ) as the objective function, the genetic algorithm (GA) is applied to inverse θ s ; if ψ(v) or ψ(Q) is used as the objective function separately, the response surface of α-n (Figure 5e), n-K s (Figure 5f), and α-K s (Figure 6d) shows that synchronous inversion of α, n, and K s is feasible, but the solutions are not unique. However, the stability of the inversion results can be improved if both ψ(v) and ψ(Q) are used as the objective functions simultaneously. The detailed objective function and solution method used in the inversion process of soil hydraulic characteristic parameters are the same as our previous paper (Yi-bo Li et al., 2018).

Inverse Solutions
A total of 12 case studies were conducted using four typical RETC soils (clay, silt, loam, and sand) with three initial water content levels, to validate the feasibility of the method. The 12 case results of the parameter inversion are summarized in Table 3, indicating that the inverted and real values were similar for each parameter.  Table 4 lists the real and inverted parameter values that were employed to draw SWRCs and SWCCs, respectively (Figure 7). SWRCs were labeled as a1, b1, c1, and d1; SWCCs were labeled as a2, b2, c2, and d2. There are small differences between the real and inverted parameter values of K s , n, and α, however, the real and inverted SWCCs and SWRCs are the same, indicating that the inverse approach is robust and effective. The error evaluation metrics NS, PBIAS, RMSE, and MAE were used to estimate the difference between the inverted and real curves, in order to further assess the robustness and accuracy of the proposed inverted approach. An average of 50 points was selected from the conductivity curves, K(θ) and suction curves, h(θ), followed by a calculation of the corresponding 50 soil water content levels. The calculated error evaluation indicators are shown in Table 4. The MAE, RMSE, and PBIAS values for the SWRCs were extremely small (0.00086-0.00164 cm 3 ·cm −3 , 0.00101-0.00192 cm 3 ·cm −3 , and 0.0102-0.3774%, respectively); the values for conductivity curves were also small (0.00501-0.01468 cm 3 ·cm −3 , 0.00800-0.02519 cm 3 ·cm −3 , and 2.1942-7.3756%, respectively), and the value of NS was relatively large (approximately 1.0). These values verify that the inverted method developed in this study is applicable to the inverse modeling of soil hydraulic characteristics [41,42]. Furthermore, the values indicate that the inverse method may be used to estimate the robustness and precision of SWRCs and SWCCs.

Inverse Solutions with Measurement Errors
There are always errors during the process of instrumentation, calibration, and other factors in practice. Hence, both the system and random errors in the initial and final water content levels, as well as the random error in the cumulative infiltration, were considered to evaluate. Table 5 shows the stability of the inverse solution. Firstly, the two system errors were set as +0.02 for θ initial and −0.02 for θ final , and the random error was set as 0.05 of the standard deviation (5.0%). The two random errors were then superimposed on the cumulative infiltration after being set as 0.02 and 0.05 of the standard deviation. All the inversion computations performed in the error-free inverse method were repeated for the four typical soils with θ initial (40% effective saturation).

Inverse Solution Analysis Based on Initial and Final Water Content
The inverse parameter values (K s , n, α, and θ s ) of the four typical soils are shown in Table 5. The system and random errors were superimposed on the error-free data, leading to only small deviations from the real parameters. It can be seen from Table 5 that most values of the inverse parameters were close to the real values, indicating that the approach is appropriate for practical application.
SWRCs and SWCCs were plotted to estimate the effective measurement error of the inverted value, as shown in Figures 8 and 9. The results of SWRCs and SWCCs show that the effective measurement error is still acceptable though there are small differences between the system and random error curves and the real value curves. Table 6 displays the error evaluation metrics NS, PBIAS, and RMSE that were used to estimate the difference between the inverted results with measurement errors and the real values. An average of 50 points was selected from the conductivity curves, K(θ) and the suction curves, h(θ), followed by a calculation of the corresponding 50 soil water content levels. Table 6 presents the results of error evaluation indicators: minimal RMSEs (0.0053-0.0346 from h(θ), 0.0004-0.0216 from K(θ)); small PBIASs (0.3181-14.5628 from h(θ) 0.1721-9.3637 from K(θ)), and large NSs (0.9322-0.9979 from h(θ), 0.9697-1.0000 from K(θ)). The results suggest that the approach is applicable to the estimation of soil hydraulic properties.

Inverted Solution Analysis Based on Cumulative Infiltration
Random errors only occurred during the process of cumulative infiltration; therefore, two standard deviations of 0.02 and 0.05 for random errors were analyzed in the inverted values (Table 7). Small deviations were found between the inverted values calculated from the values with measurement errors and the real values. Figure 10 shows SWRCs and SWCCs based on the inverse values with two random errors. The inverse method is robust and reasonable for the inversion of soil hydraulic parameters, due to that, there are small differences between the inverted and real values of these curves. Table 8 details the assessment criteria NS, PBIAS, and RMSE for the inverted and real values. The robustness and effectiveness of the inverse method are represented by large NS, small PBIAS, and RMSE.

Inverse Solution Based on Experimental Infiltration Data
Comparisons of parameter estimation results for the sandy loam are depicted in Table 9. Two h(θ) curves and K(θ) curves with hydraulic parameters inversed from the laboratory experiments are presented in Figure 11a1,a2,b1,b2, respectively. Figure 11a,b represent the results of two bulk densities (1.30 and 1.35 g cm −3 ) with three initial soil water contents. The three soft dotted lines (circle-18.92%, triangle-26.17%, and square-31.74%) were determined separately by using the VEGA-PSO inverse method, and the solid line was then determined from all three datasets from the laboratory experiments and using the same inverse method. Figure 11 illustrates that despite the small differences between the three sets of parameter values and the solid lines, they are close to each other. Table 10 displays RMSE, PBIAS, and NS for each inverse value and unified value. The robustness and effectiveness of the inverse method are represented by large NS, small PBIAS, and RMSE. Therefore, the results proved that the method proposed in this paper can be applied in laboratory experiments to invert the soil hydraulic parameters. Meanwhile, it is more intuitive to compare the two curves, h(θ) curves and K(θ) curve, to check the accuracy of the inversion than to compare the values of the individual soil hydraulic parameters individually.

Comprehensive Analysis of the Inverse Solution of Experimental Results
To further validate the use of the inverse solution with laboratory experimental data, the cumulative infiltration and final water content were compared between the inverse values (unified values) and the measured values. Figure 12 depicts comparisons of the cumulative infiltration for bulk densities of 1.30 and 1.35 g cm −3 . Small deviations can be seen among the six graphs, which may be due to errors in the experiments. Figure 13 shows the distribution of final water content in the six experiments.             differences between the three sets of parameter values and the solid lines, they are close to each other. Table 10 displays RMSE, PBIAS, and NS for each inverse value and unified value. The robustness and effectiveness of the inverse method are represented by large NS, small PBIAS, and RMSE. Therefore, the results proved that the method proposed in this paper can be applied in laboratory experiments to invert the soil hydraulic parameters. Meanwhile, it is more intuitive to compare the two curves, h(θ) curves and K(θ) curve, to check the accuracy of the inversion than to compare the values of the individual soil hydraulic parameters individually.
(b1) h(θ) 1.35 g cm −3 ( b2) K(θ) 1.35 g cm −3 Figure 11. Comparison of soil water conductivity curves and SWRCs of the three experimental solutions. Figure 11. Comparison of soil water conductivity curves and SWRCs of the three experimental solutions.  To further validate the use of the inverse solution with laboratory experimental data, the cumulative infiltration and final water content were compared between the inverse values (unified values) and the measured values. Figure 12 depicts comparisons of the cumulative infiltration for bulk densities of 1.30 and 1.35 g cm −3 . Small deviations can be seen among the six graphs, which may be due to errors in the experiments. Figure 13 shows the distribution of final water content in the six experiments.
(a1) 18 Figure 12. Comparison of cumulative infiltration for three initial water content levels. Figure 12. Comparison of cumulative infiltration for three initial water content levels.

Conclusions
This study investigated the inverse determination of soil hydraulic properties based on the van Genuchten-Mualem model and data from simultaneously measured cumulative infiltration and final water content. The proposed two-step inverse estimation method Figure 13. Comparison of final water distribution for three initial water content levels.

Conclusions
This study investigated the inverse determination of soil hydraulic properties based on the van Genuchten-Mualem model and data from simultaneously measured cumulative infiltration and final water content. The proposed two-step inverse estimation method was employed to compare numerically simulated values and measurements from laboratory experiments, indicating that the data were consistent. Moreover, new formulations of α/K s and v/K s were applied in the inverse procedure to reduce the scope of K s , increasing the efficiency of the method. After the inverse procedure was performed, the GA with θ final was used to determine the saturated water content θ s . Finally, the soil characteristic parameters K s , n, and α were accurately inversed by using the VEGA-PSO-based multiobjective optimization method according to the infiltration rate, final water content, and cumulative infiltration.
In terms of validating the accuracy of the model, no separate comparison of the values of each soil hydraulic parameter was made, while considering that the parameters previously influenced each other, and two curves h(θ) and K(θ) were plotted for accuracy analysis. Additionally, the model was verified in the case of measurement errors (system error and random error) that exist in the practical situation.
The general agreement between the numerical simulations and inverse solutions for h(θ) and K(θ) suggested that the proposed inverse method is very effective. The results obtained directly from laboratory experiments were also found to agree with the inverse solutions. The robustness and practicability of the proposed method were validated given that the accuracy of results was acceptable. However, further investigation is required for the inhomogeneity of soil water content distribution and soil texture. Moreover, the estimation of soil hydraulic parameters on field plot and large spatial scale needs more thorough research. Li, Y.-B.; Liu, Y.; Ma, X.-Y.