Stable expression and control parameters in the numerical simulation of unsaturated flow

The Richards’ equation describes the flow phenomenon in unsaturated porous media and is essential to hydrology and environmental science. This study evaluated the numerical stability of two different forms of the Richards’ equation. Sensitivity analyses were performed to investigate the control parameters of the equation. The results show that the h-form Richards’ equation has better applicability for calculating variable saturation flows than the θ-form Richards’ equation. For the h-form Richards’ equation, the hydraulic conductivity of the soil in the low-suction range and the specific moisture capacity in the high-suction range primarily influenced the solution. In addition, sensitivity analyses indicated that the saturated hydraulic conductivity, initial condition, and air-entry pressure have a higher sensitivity to the simulation results than the saturated water content, rainfall intensity, and decline rate of hydraulic conductivity. Moreover, their correctness needs to be guaranteed first in numerical simulations. The research findings can provide a helpful reference for improving the reliability of numerical simulations of unsaturated flows.


Introduction
Soil moisture prediction is significant to hydrology, environmental science, agriculture, and pollutant treatment [1]. With the development of computer technology, numerical simulations have gradually become the primary method for studying soil water movement. The Richards' equation [2] describes the flow phenomenon in unsaturated porous media under the action of gravity and capillary forces. It has been widely used in numerical models of unsaturated flows and is essential to hydrology and environmental science [3,4].
The solution of the Richards' equation requires determining soil water parameters, such as hydraulic diffusivity D, hydraulic conductivity k, and specific moisture capacity C, which comprehensively reflect the water holding and transport capacity of the soil [5]. For decades, considerable effort has been devoted to determining or inferring soil water parameters [6][7][8][9][10][11]. The hydraulic diffusivity D of soil is the ratio of the hydraulic conductivity k to the specific moisture capacity C, reflecting the soil porosity, pore size distribution, and water conductivity, and affects water movement in the soil. The horizontal soil column method proposed by Bruce and Klute [12] is the most used method for measuring D. The limitation of this method is that the hydraulic diffusivity obtained typically exhibits significant discreteness. In particular, at a high water content, it is difficult to measure the hydraulic diffusivity accurately [13,14]. The curve of the specific moisture capacity C derived from the conventional soil water characteristic curve models in geotechnical engineering (e.g., Brooks and Corey model [15], van Genuchten model [16], and Fredlund and Xing model [17]) is typically bell-shaped. C is either small or zero when the soil is near dry and wet. Because the hydraulic diffusivity D = k·dψ/dθ, it is close to infinity when the soil is saturated; that is, when the θform Richards' equation is used for the numerical calculation, D(θ s ) (i.e., D at saturation) cannot be determined, making the θ-form Richards' equation unable to simulate the partially saturated soil water flow. The conventional measure to solve this problem is to artificially specify the D(θ s ) value to avoid infinite hydraulic diffusivity [3,[18][19][20][21]. Most studies on numerical solutions of the Richards' equation are based on various post-processing [22][23][24][25][26][27][28], and the influence of model parameters on the numerical results is still unclear.
Parameter measurement errors are the primary sources of uncertainty in numerical models, significantly affecting the reliability of soil moisture predictions [29]. Sensitivity analyses have been used to study the influence of parameter changes in the model on the output evaluate the importance of each parameter to the simulation results, improving the reliability [30][31][32]. There are many types of sensitivity analysis methods. Among them, the Morris screening method [33] is an effective global sensitivity analysis method and is favored by researchers owing to its small computational load [34][35][36]. King and Perera [37] used the Morris method to study an urban water supply system and analyzed the importance of the input variables. Ren et al. [38] used the Morris method to study the primary factors affecting the embankment temperature field using model parameters based on a hydrothermally coupled model. Yi et al. [39] used the Morris method to study the sensitivity of the sample size and parameter perturbation range to different output indicators of the model. Moreover, they analyzed the primary factors controlling nutrient migration in a complex three-dimensional water quality model. In conclusion, owing to the errors and uncertainties in the input parameters, it is necessary to conduct a sensitivity analysis as an essential prerequisite in the numerical simulation of moisture migration. However, few studies have analyzed the primary controlling factors in the numerical models of water migration.
In this study, a numerical method was used to analyze the influence of hydraulic diffusivity on the seepage process when describing soil water infiltration using two forms of the Richards' equation. The numerical stability of the two forms of Richards' equation was evaluated, and suggestions were provided for their use. In addition, taking the pressure head-form Richards' equation as an example, the Morris sensitivity analysis method was used to study the sensitivity of the six parameters in the numerical model, providing a useful reference for improving the reliability of the numerical simulation of water migration. The Richards' equation is typically used to describe unsaturated flow phenomena in porous media (e.g., soil). The original form of the Richards' equation, referred to as the mixed Richards' equation (or mixed form), contains two variables (i.e., pressure head h and water content θ). Although the Richards' equation is easy to deduce, it is a challenging equation to solve accurately in hydrology owing to its high nonlinearity [3]. Converting the mixed Richards' equation into a form with only one unknown simplifies the problem, whether it is to obtain an analytical solution or a numerical solution. Therefore, the single-variable Richards' equation (the pressure head form or water content form) is typically used to simulate water movement in the saturated-unsaturated zone [40].
The Richards' equation describes the one-dimensional horizontal infiltration problem of a soil column. The Richards' equation in the form of a pressure head (hform) and the corresponding initial and boundary conditions are as follows: where h is the soil suction, C(h) is the specific moisture capacity of the soil, k(h) is the hydraulic conductivity, h 0 is the initial suction of the soil, and x and t are the horizontal coordinates and time, respectively. The horizontal coordinates take the left boundary as the origin and the horizontal right as the positive direction. The Richards' equation in the form of water content (θform) and the corresponding initial and boundary conditions for one-dimensional horizontal flows are as follows: (2) where θ is the volumetric water content, D(θ) is the ratio of the hydraulic conductivity to the specific moisture capacity, called the hydraulic diffusivity for unsaturated soil, θ 0 is the initial soil water content, and θ s is the saturated water content of soil (i.e., the water content of soil at the boundary).

Method for solving the Richards' equation
COMSOL Multiphysics is a powerful software with strong multiphysical coupling and nonlinear differential equation solving capabilities. In this study, through the secondary development of the convection-diffusion equation module in COMSOL, numerical models of saturated-unsaturated water movement based on the θform and h-form Richards' equations were established. The expression of the convection-diffusion equation and boundary condition provided in COMSOL is as follows: ∇ where d a is the damping or mass coefficient, u is a variable, is a differential operator, [∂/∂x] is a onedimensional problem, c is the diffusivity, β is the convection coefficient, f is the source term, r is the boundary condition parameter, Ω is the calculation domain, and Γ is the boundary of the calculation domain.
∂θ ∂t 3 Numerical stability analyses for unsaturated flows

Description of numerical model and plan
In the numerical tests, a one-dimensional soil column with a length of 0.15 m was used, and the influence of gravity was ignored. The soil column was horizontal, as shown in Fig. 1. To avoid numerical oscillations caused by the strong nonlinearity of the hydraulic characteristic function, the grid was divided into linear elements of 0.001 m. The boundary condition of the model was set at the left boundary. The boundary conditions remained unchanged during the water infiltration from left to right. The soil column was uniform with the same initial conditions (suction ψ and water content θ) and hydraulic characteristics. Table 1 lists the numerical test plans. The soil hydraulic characteristics used in the numerical tests were obtained from Fredlund and Xing [17]. The soil water characteristic curve and soil hydraulic conductivity function were also used, described by the van Genuchten model [16] and Mualem model [41], respectively, as shown in Figs. 2(a) and 2(b). The Dirichlet boundary conditions were used in the numerical tests: (1) in the h-form Richards' equation, u = 0 kPa; (2) in the θ-form Richards' equation, u = 0.322 (volumetric water content), which does not cause ponding on the surface of the soil column. In this study, the initial soil suction was 27000 kPa (or the initial water content was 0.110), and the soil was relatively dry. The hydraulic diffusivity D can be changed by adjusting the hydraulic conductivity k or specific moisture capacity C. The influence of these two variables on the seepage calculation was analyzed. When one variable is changed, the other variable is fixed: (1) change the hydraulic conductivity k and keep the specific moisture capacity C unchanged; and (2) change the specific moisture capacity C and keep the hydraulic conductivity k unchanged. To analyze the influence of hydraulic diffusivity on the seepage process, the hydraulic diffusivity is treated as follows, as shown in Figs. 2(c) and 2(d): • Both T1 (ABCD) and H1 (ABCD) represent experimental data. T1 is the relationship between volumetric water content and hydraulic diffusivity, and H1 is the relationship between suction and hydraulic diffusivity.
• H2 and H6 (FBCD) had a 10-fold reduction in hydraulic diffusivity at ψ = 0.1 kPa (corresponding to θ = 0.322 in T1) compared with H1. The hydraulic diffusivity was adjusted by changing the specific moisture capacity in H2 and hydraulic conductivity in H6.
Notes: a) D s (m 2 /s) represents the saturated hydraulic diffusivity of the soil; b) D 0 (m 2 /s) represents the initial hydraulic diffusivity; c) variable is the hydraulic characteristic function changed in different numerical tests, which is C or k in the h-form Richards' equation and D in the θ-form Richards' equation; d) ψ 0 is the initial soil suction and θ 0 is the initial water content; the hydraulic characteristic functions used in e) H1 and f) T1 were the experimental data from Fredlund and Xing [17].
• H3 and H7 (EBCD) had a 10-fold increase in hydraulic diffusivity at ψ = 0.1 kPa (corresponding to θ = 0.322 in T1) compared with H1. The hydraulic diffusivity was adjusted by changing the specific moisture capacity at H3 and hydraulic conductivity at H7.
• H4 and H8 (ABCH) had a 10-fold reduction in hydraulic diffusivity at ψ = 27000 kPa (corresponding to θ = 0.110 in T1) compared with H1. The hydraulic diffusivity was adjusted by changing the specific moisture capacity in H4 and hydraulic conductivity in H8.
• H5 and H9 (ABCG) had a 10-fold increase in hydraulic diffusivity at ψ = 27000 kPa (corresponding to θ = 0.110 in T1) compared with H1. The hydraulic diffusivity was adjusted by changing the specific moisture capacity in H5 and hydraulic conductivity in H9.

θ-form Richards' equation
Figure 2(c) shows the variation in the hydraulic diffusivity in the numerical tests T1, T2, and T3 when the soil was close to saturation. The value of the hydraulic diffusivity changed significantly between volumetric water contents of 0.313 and 0.322. Figure 3(a) shows the water content profile of the soil column at t = 1000 s. In the numerical test, the soil moisture infiltration distance was determined by the position of the wetting front, and the characteristic water content θ d = 0.12 was used, i.e., when the volumetric water content at a certain position in the soil column increased to 0.12, it was determined that the wetting front reached this position. Figure 3(b) shows the relationship between the advancing distance and time of the wetting front in the three numerical tests. The advancing distance of the wetting front increases with increasing hydraulic diffusivity in the saturated section. At t = 300 s, the advancing distances of the wetting front in the numerical tests of T1, T2, and T3 were 0.039, 0.031, and 0.056 m, respectively. The advancing distances of the wetting front in numerical tests T2 and T3 were 0.79 times and 1.44 times that of T1, respectively. Figure 2(c) also shows the variation in the hydraulic diffusivity in numerical tests T4 and T5 at the beginning of soil infiltration (i.e., high-suction range). The hydraulic diffusivity varied between volumetric water contents of 0.070 and 0.150. Figure 4 shows the water content profiles of the soil columns in numerical tests T1, T4, and T5 at t = 1000 s. The change in hydraulic diffusivity of the residual section did not affect the water seepage process.

h-form Richards' equation
Figure 2(d) shows the change in hydraulic diffusivity in numerical tests H2, H3, H6, and H7 when the soil was near saturation. The hydraulic diffusivity changed significantly between 0.1 and 80 kPa soil suction. Figure 5 shows the specific moisture capacity data used in the numerical tests H1, H2, and H3 when the hydraulic diffusivity is changed by the specific moisture capacity. Figure 6 shows the hydraulic conductivity function used in the numerical tests H1, H6, and H7 when the hydraulic diffusivity is changed by hydraulic conductivity. Figure 7(a) shows the suction profile of the soil column at t = 1000 s. Figure 7(b) shows the time history of the advancing distance of the wetting front when the characteristic suction ψ d = 25000 kPa is selected. As shown in Fig. 7, for soils that are close to saturation, changing the hydraulic diffusivity using different methods will have different influences on the seepage process. The change in the hydraulic diffusivity did not affect the seepage process when the specific moisture capacity was used as the control variable for the calculation. As shown in Fig. 7(b), the advancing distances of the wetting front in numerical tests H1, H2, and H3 were the same. However, hydraulic diffusivity changed the seepage process when hydraulic conductivity was used as the control variable. At 300 s, the advancing distances of the wetting front in numerical tests H1, H6, and H7 were 0.038, 0.032, and 0.054 m, respectively. The advancing distances of the wetting front in numerical tests H6 and H7 were 0.84 times and 1.42 times that of H1, respectively.   Figure 2(d) shows the change in hydraulic diffusivity in numerical tests H4, H5, H8, and H9 at the beginning of soil infiltration (i.e., the high-suction range). The hydraulic diffusivity changed significantly between 7000 and 30000 kPa. Figure 8 shows the specific moisture capacity data used in the numerical tests H4 and H5 when the hydraulic diffusivity is changed by the specific moisture capacity. Figure 9 shows the hydraulic conductivity functions used in numerical tests H8 and H9 when the hydraulic diffusivity is changed by the hydraulic conductivity. Figure 10(a) shows the suction profile of the soil column at t = 1000 s. Figure 10(b) shows the time history of the advancing distance of the wetting front when the characteristic suction ψ d = 25000 kPa is selected. As shown in Fig. 10, at the beginning of soil infiltration, the change in hydraulic diffusivity only had a slight impact on the seepage process when hydraulic conductivity was used as the control variable. At t = 300 s, the advancing distances of the wetting front in numerical tests H8 and H9 were 0.037 and 0.040 m, respectively. The advancing distances of numerical tests H8 and H9 were 0.98 and 1.05 times that of H1, respectively. When the specific moisture capacity was used as the control variable, the hydraulic diffusivity changed the seepage process of soil water. At t = 300 s, the advancing distances of the wetting front in numerical tests H4 and H5 were 0.023 and 0.044 m, respectively. The advancing distances of numerical tests H4 and H5 were 0.61 and 1.16 times that of H1, respectively.

Comparison between two expressions
The h-form and θ-form Richards' equation were numerically calculated to analyze the influence of hydraulic diffusivity on soil water infiltration.
For the θ-form Richards' equation, when the soil was close to saturation, the advancing speed of the wetting front increased with increasing D. When the hydraulic diffusivity D increased by 10 times, the advancing distance of the wetting front was 1.44 times that before the change. When the hydraulic diffusivity D decreased to one-tenth, the advancing distance of the wetting front was 0.79 times that before the change. Changes in hydraulic diffusivity at the beginning of soil infiltration (high-suction range) had negligible effects on the flow process.
For the h-form Richards' equation, different treatment methods had different effects on the seepage process   when hydraulic diffusivity was changed. When the soil was near saturation, the change in the specific moisture capacity C did not affect the seepage process. Despite a change in the specific moisture capacity C by 10 times when the soil was saturated (ψ = 0.1 kPa), the water seepage process remained unchanged. However, when the soil was near saturation, the hydraulic conductivity k changed the seepage process, and the advancing speed of the wetting front increased with increasing k. When the hydraulic conductivity k increased by 10 times, the advancing distance of the wetting front was 1.42 times that before the change. When the hydraulic conductivity k decreased by 10 times, the advancing distance of the wetting front was 0.84 times that before the change. However, at the beginning of the soil infiltration (highsuction range), the change in hydraulic conductivity k had a negligible impact on the seepage process. The advancing speed of the wetting front decreased with an increase in the specific moisture capacity C. When the specific moisture capacity C was 10 times larger, the advancing distance of the wetting front was 0.61 times before the change. When the specific moisture capacity C was 10 times smaller, the advancing distance of the wetting front was 1.16 times before the change.
In summary, when using the θ-form Richards' equation, the correctness of the hydraulic diffusivity data in a lowsuction range (i.e., the soil is close to saturation) should be ensured. By contrast, when using the h-form Richards' equation, it is necessary to control the measurement error of the hydraulic conductivity function in the low-suction range and specific moisture capacity data in the highsuction range.

Comparison of numerical results with existing benchmarks
This section compares the numerical results with existing benchmarks to verify the accuracy of the numerical method. Crank [42] deduced the following infiltration formula for soil: where Q is the cumulative infiltration, is the weighted mean diffusivity, and α is the soil parameter, set to 25.16 in this calculation. The soil used in tests T1-T5 was reanalyzed using Eq. (6). Table 2 lists the infiltration results after 300 s. The results calculated using Eq. (6) are consistent with the numerical results of the θ-form Richards' equation. The infiltration process is susceptible to diffusivity in the low-suction section. This is because in the low-suction section, the water content θ and diffusivity D increase rapidly (as shown in Fig. 2(c)), increasing the weighted mean diffusivity and cumulative infiltration Q; that is, D(θ s ) is crucial for the numerical calculation of the θ-form Richards' equation. However, D(θ s ) cannot be determined when the soil reaches saturation during the wetting process, resulting in different solutions to the θ-form Richards' equation in the numerical calculation (i.e., failure to solve the equation), as shown in Fig. 11(a). However, whether k or C causes a change in the infiltration cannot be verified using Eq. (6), because the soils used in tests H3 and H7 had the same. Tracy [43] provided an analytical solution for transient unsaturated seepage to test the accuracy of the finite element method.
where h is the infiltration distance, α and c are soil parameters, h r is the initial suction of the soil, L is the height of the soil column, λ k = kπ/L, and µ k = (α 2 /4+λ k 2 )/c. The soil used in tests H1, H3, and H7 was reanalyzed using Eq. (7) to calculate the infiltration distance. Figure 11(b) compares the distance of the wetting front obtained by the numerical solution of the hform Richards' equation and analytical solution [43]. The numerical solution is consistent with the analytical solution. The advancing distance curves of H1 and H3 coincided, whereas the infiltration speed of H7 was higher. The hydraulic conductivity k in the low-suction range was the primary factor controlling the seepage. The existing measurement method of k in the low-suction range is established [44][45][46]; therefore, the h-form Richards' equation is suitable for simulating the partial saturation movement.

Method of sensitivity analyses
Sensitivity analysis is a method for studying the allocation of the uncertainty of the model output to the uncertainty of the model input factors [47]. Sensitivity analysis of the parameters in the model is an essential step when solving seepage problems using numerical methods. It effectively evaluates the contribution rate and degree of influence of the model parameters on the seepage process, improving the reliability of the numerical model. The Morris screening is a widely used sensitivity analysis method [48]. In this method, only one parameter x i in the model is changed during the analysis, with the remaining parameters fixed. The model outputs the objective function y(x) = y(x 1 ,x 2 ,x 3 ,...,x n ) using the indicator e i to determine the influence of the parameter change on the model output. The formula for e i is as follows: where y and y i are the output values of the model before and after the parameter change, and x and x i are the parameters before and after the change, respectively. In this study, six parameters in the numerical model were selected for sensitivity analyses: saturated water content θ s , saturated hydraulic conductivity k s , initial condition ψ 0 , air-entry value ψ b , rainfall intensity v, and hydraulic conductivity decline rate. Owing to the different units of each parameter and significant difference in their value range, three parameters (i.e., saturated hydraulic conductivity k s , initial condition ψ 0 , and airentry value ψ b ) were standardized using Eq. (9) to facilitate comparison. The other three parameters (i.e., saturated water content θ s , rainfall intensity v, and hydraulic conductivity decline rate) were standardized using Eq. (10).
Francos et al. [49] proved that when analyzing the sensitivity of the parameters, the calculation accuracy of Eq. (8) is higher if the parameters change at a fixed step. The sensitivity coefficient S i is the average value of the calculation results after multiple parameter disturbances. The modified Morris screening method is as follows: /(n − 1), (11) where Y i is the output value of the ith operation of the model, Y i+1 is the output value of the (i+1)th operation of the model, Y 0 is the output value of the model when the input parameter is the benchmark parameter, P i is the percentage change in the parameter value relative to the benchmark parameter during the ith operation of the model, P i+1 is the percentage change in the parameter value relative to the benchmark parameter during the (i+1)th operation of the model, and n is the number of model operations.

Description of numerical model
The numerical tests in this section used the same numerical model as that in Section 3. The model length was set to 25 m to fully observe the advancing distance of the wetting front. Under rainfall conditions, the infiltration boundary can be switched between the Neumann boundary (i.e., velocity boundary) and Dirichlet boundary (i.e., pressure boundary) according to the topsoil conditions [50]. In COMSOL, a complementary smoothing function can define the mixed boundary condition, expressed as where n is the outer normal vector of the boundary, k s is the saturated hydraulic conductivity, ρ is the density of water (1000 kg/m 3 ), p is the pore water pressure, z is the ordinate, g is the acceleration of gravity, and v is the rainfall intensity. H b = z + h b is the external total water head, where h b is the water depth, takin as 0.01 m in the calculation. H = z + p/(ρg) is the total water head at the boundary, R b is the external resistance (R b = 1000k s ), and α and β are complementary smooth functions, as shown in Eq. (13): In the numerical calculation, the conversion between the two infiltration boundary conditions was controlled according to the pore water pressure p of the surface soil. The infiltration boundary condition was set as the velocity boundary when the pore water pressure p was less than 0. When the pore water pressure p was greater than or equal to 0, the infiltration boundary condition was set as the pressure boundary.

Scheme of numerical tests
The unsaturated seepage process was described through numerical tests using the h-form Richards' equation. The sensitivity of various parameters in the seepage model was analyzed. According to the classification standard of rainfall grade in meteorology, two rainfall modes (longterm light rain and short-term heavy rain) were selected. Three typical soils (clay, silt, and sand) were analyzed in the test. Table 3 shows the model parameters in the test, in which the rainfall intensity was disturbed up and down in steps of 20%. Table 4 shows the ranges of the soil parameters [14] and their values from the numerical tests. Each parameter was disturbed up and down in steps of 20%. Figure 12 shows the benchmark values of the soil moisture characteristic function used in the numerical tests.

Results and analyses
The numerical test considered the advancing distance of the wetting front as the model output. Figures 13 and 14 show the test results. Within the range of parameter changes, the three parameters of rainfall intensity v, airentry pressure ψ b , and saturated hydraulic conductivity k s are positively correlated with the advancing distance of the wetting front; that is, the advancing distance of the wetting front increases with an increase in the parameters. By contrast, the saturated water content θ s , hydraulic conductivity decline rate, and initial condition ψ 0 were negatively correlated with the advancing distance of the wetting front; that is, the advancing distance of the wetting front decreased with an increase in the parameters.
The sensitivity of each parameter was calculated using Eq. (9). Figure 15 shows the calculation results. For the three typical soils under the two rainfall modes, the saturated hydraulic conductivity k s , initial condition ψ 0 , and air-entry pressure ψ b are highly sensitive. The initial condition ψ 0 of clay is the most sensitive parameter under the long-term light rain mode. The initial condition ψ 0 of for sand and silt was the most sensitive parameter under both rainfall modes. Therefore, it is essential to determine the initial field when analyzing the water migration caused by rainfall. Under the short-term heavy rain mode, the sensitivity of initial condition ψ 0 of clay and silt was significantly reduced. This is because when the rainfall intensity v is greater than the saturated hydraulic conductivity k s of soil (the benchmark saturated hydraulic conductivities of clay and silt are 3e−9 and 3e−7 m/s, respectively, and the rainfall intensity is 2e−6 m/s), the   degree of saturation of the soil increases rapidly in a short time. Consequently, the influence of the initial condition ψ 0 on the seepage process is reduced. The saturated water content θ s , rainfall intensity, and decline rate of hydraulic conductivity were less sensitive than the other parameters. Therefore, the accuracy of the saturated hydraulic conductivity k s , initial condition ψ 0 , and air-entry value ψ b should be ensured first in the numerical calculation.

Engineering case analysis
In this section, a two-dimensional reservoir slope is used to verify the validity of the sensitivity analysis results presented in Subsection 4.2. Figure 16 shows the numerical model of the engineering slope and typical section. Most of the slope was covered by silt. The natural slope was generally between 12° and 28°. The slope was 16 m high and 50 m wide. The bottom of the slope was a bedrock. The groundwater level was 6 m above the bedrock. A geometric model with a length of 54 m and height of 28 m was established to reduce the influence of the boundary [20]. The two sides of the model were symmetrical boundaries, and the bottom side was an impervious boundary simulating the bedrock. The standardized model parameters were increased by 20% to analyze the influence of the model parameters on the calculation results. Table 5 lists the initial and adjusted parameters of the model.
According to the climatic conditions, a rainfall intensity of 38 mm/h and rainfall duration of 10 h were used to simulate typical local rainfall conditions. Figure 17(a) shows the saturation distribution of the slope after rainfall, simulated based on geological exploration data. The rainfall infiltration depth in the section 10 m from the left boundary was analyzed (as shown by the dotted line in Fig. 17(a)). Figure 17(b) illustrates the absolute difference between the infiltration depth after the parameter change and the initial infiltration depth. The saturated hydraulic conductivity k s , initial condition ψ 0 , and air-entry pressure ψ b are the primary factors affecting the results. By contrast, the other parameters that have little impact on the model output can be simplified under appropriate assumptions in engineering applications. Figure 17(c) shows the distribution of the slope saturation after the change in each parameter. The engineering slope analysis results verify the validity of the conclusions presented in Subsection 4.2.

Conclusions
In this study, COMSOL software was used to simulate the process of unsaturated flow based on the h-form and θ-form Richards' equations, and the influence of hydraulic diffusivity on the seepage process was studied. In addition, sensitivity analyses of the parameters in the     numerical model were performed using the h-form Richards' equation. The main conclusions are as follows. 1) For the θ-form Richards' equation, the advancing speed of the wetting front increased with increasing saturated hydraulic diffusivity D s . By contrast, the changes in hydraulic diffusivity D in the high-suction range (i.e., at the beginning of infiltration) did not affect the seepage process. For the h-form Richards' equation, the advancing speed of the wetting front increased with increasing saturated hydraulic conductivity k s . The specific moisture capacity C in the high-suction range (i.e., at the beginning of infiltration) also affected the seepage process. Therefore, when using the θ-form Richards' equation to calculate the seepage process, the calculation results were susceptible to the hydraulic diffusivity in the low-suction range of the soil, and there was significant numerical uncertainty in the calculation. Therefore, the θ-form Richards' equation was not recommended for seepage calculations when soil was close to saturation. The h-form Richards' equation had good applicability in the calculation of variable saturation flow. The seepage process was primarily controlled by the hydraulic conductivity in the low-suction range and specific moisture capacity in the high-suction range.
2) In the numerical model, the different parameters had different effects on the infiltration process of water in the soil. The three parameters of rainfall intensity v, air-entry pressure ψ b , and saturated hydraulic conductivity k s were positively correlated with the seepage process. The advancing speed of the wetting front increased with an increase in these parameters. However, the saturated water content θ s , hydraulic conductivity decline rate , and initial condition ψ 0 negatively affected the water seepage process, i.e., the advancing speed of the wetting front decreased with an increase in these parameters. These findings provided valuable information for the revision of moisture-transport models.
3) Sensitivity analyses of the six parameters indicated that under different rainfall modes (e.g., long-term light rain and short-term heavy rain) the soil parameters of saturated hydraulic conductivity k s , initial condition ψ 0 , and air-entry pressure ψ b were highly sensitive to the simulation results. Therefore, their correctness needed to be guaranteed first in numerical simulations. By contrast, other parameters (i.e., saturated water content θ s , rainfall intensity v, and decline rate of hydraulic conductivity) with relatively low sensitivity could be simplified under appropriate assumptions. The selection of high-sensitivity parameters provided a basis for determining the priority of parameters in numerical simulations of unsaturated flows. Moreover, it provided a useful reference for improving the reliability of the numerical model of moisture migration in engineering practice.