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

Abstract

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.

Graphical abstract

Keywords

Richards’ equation / sensitivity analysis / unsaturated soil / hydraulic diffusivity / seepage simulation

1 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–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–21]. Most studies on numerical solutions of the Richards’ equation are based on various post-processing [22–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–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–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.

2 Unsaturated flows in porous media

2.1 Differential expression of unsaturated flow

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 (h-form) and the corresponding initial and boundary conditions are as follows:

(1)\{\begin{array}{rl}& C(h)\frac{\mathrm{\partial }h}{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(k(h)\frac{\mathrm{\partial }h}{\mathrm{\partial }x}\right),\\ & h=h_0,x>0,t=0,\\ & h=0,x=0,t>0,\end{array}

where h is the soil suction, C(h) is the specific moisture capacity of the soil, k(h) is the hydraulic conductivity, h₀ 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)\{\begin{array}{rl}& \frac{\mathrm{\partial }\theta }{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(D(\theta )\frac{\mathrm{\partial }\theta }{\mathrm{\partial }x}\right),\\ & \theta ={\theta }_0,x>0,t=0,\\ & \theta ={\theta }_{\text{s}},x=0,t>0,\end{array}

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, θ₀ 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).

2.2 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:

(3)\{\begin{array}{c}{d}_{\text{a}}{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}}+\mathrm{\nabla }\cdot (-c\mathrm{\nabla }u)+\beta \cdot \mathrm{\nabla }u=f,\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ u=r,\mathrm{o}\mathrm{n}\mathrm{\Gamma },\end{array}

where d_a is the damping or mass coefficient, u is a variable, \mathrm{\nabla } is a differential operator, [∂/∂x] is a one-dimensional 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.

Convert Eqs. (1) and (2) into the convection–diffusion equation provided by COMSOL.

(4)C(h)\frac{\mathrm{\partial }h}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (-k(h)\mathrm{\nabla }h)=0,

(5)\frac{\mathrm{\partial }\theta }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (-D(\theta )\mathrm{\nabla }\theta )=0.

3 Numerical stability analyses for unsaturated flows

3.1 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.

Tab.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 Fig.2(a) and Fig.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 Fig.2(c) and Fig.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.

● T2 (FBCD) had a 10-fold reduction in hydraulic diffusivity at θ = 0.322 compared with T1.

● T3 (EBCD) had a 10-fold increase in hydraulic diffusivity at θ = 0.322 compared with T1.

● T4 (ABCH) had a 10-fold reduction in hydraulic diffusivity at θ = 0.110 compared with T1.

● T5 (ABCG) had a 10-fold increase in hydraulic diffusivity at θ = 0.110 compared with T1.

● 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.

● 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.

3.2 Numerical stability analyses

3.2.1 θ-form Richards’ equation

Fig.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. Fig.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. Fig.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.

Fig.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. Fig.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.

3.2.2 h-form Richards’ equation

Fig.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. Fig.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. Fig.6 shows the hydraulic conductivity function used in the numerical tests H1, H6, and H7 when the hydraulic diffusivity is changed by hydraulic conductivity.

Fig.7(a) shows the suction profile of the soil column at t = 1000 s. Fig.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.

Fig.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. Fig.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. Fig.9 shows the hydraulic conductivity functions used in numerical tests H8 and H9 when the hydraulic diffusivity is changed by the hydraulic conductivity.

Fig.10(a) shows the suction profile of the soil column at t = 1000 s. Fig.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.

3.2.3 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 (high-suction 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 low-suction 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 high-suction range.

3.3 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:

(6)\{\begin{array}{rl}& Q=2({\theta }_{\mathrm{s}}-{\theta }_0)\sqrt{\overline{D}t/\text{π}},\\ & \overline{D}=\alpha \underset{{\theta }_0}{\overset{{\theta }_{\mathrm{s}}}{\int }}{(\theta -{\theta }_0)}^{2/3}D(\theta )\text{d}\theta ,\end{array}

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). Tab.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.

(7)\{\begin{array}{rl}& h=\frac{1}{\alpha }\mathrm{l}\mathrm{n}\left({\text{e}}^{\alpha h_{\text{r}}}+h_{\text{t}}(\mathit{\text{z}},t)+h_{\text{s}}(\mathit{\text{z}})\right),\\ & h_{\text{s}}(\mathit{\text{z}})=\left(1-{\text{e}}^{\alpha h_{\text{r}}}\right)\left(1-{\text{e}}^{-\alpha \mathit{\text{z}}}\right)/\left(1-{\text{e}}^{-\alpha L}\right),\\ & h_{\text{t}}(\mathit{\text{z}},t)=\frac{2\left(1-{\text{e}}^{\alpha h_{\text{r}}}\right)}{Lc}{\text{e}}^{\alpha (L-\mathit{\text{z}})/2}\\ & \cdot \sum _{k=1}^{\mathrm{\infty }}{(-1)}^k\left(\frac{{\lambda }_k}{{\mu }_k}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\lambda }_k\mathit{\text{z}}\right)\text{e}-{\mu }_{k}t,\end{array}

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 = (α²/4+λ_k²)/c. The soil used in tests H1, H3, and H7 was reanalyzed using Eq. (7) to calculate the infiltration distance. Fig.11(b) compares the distance of the wetting front obtained by the numerical solution of the h-form 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–46]; therefore, the h-form Richards’ equation is suitable for simulating the partial saturation movement.

4 Primary factors controlling the h-form Richards’ equation

4.1 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₁,x₂,x₃,...,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:

(8)e_i=\frac{y_i-y}{(x_i-x)/x},

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 ψ₀, 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 ψ₀, and air-entry 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).

(9)X_i=\frac{\lg x_i-\lg x_{min}}{\lg x_{max}-\lg x_{min}},

(10)X_i=\frac{x_i-x_{min}}{x_{max}-x_{min}}.

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:

(11)S_i=\sum _{i=1}^n\frac{(Y_{i+1}-Y_i)/Y_0}{(P_{i+1}-P_i)/100} /(n-1 ),

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₀ 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.

4.2 Numerical test

4.2.1 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

(12)n\times \frac{k_{\mathrm{s}}}{\rho g}\mathrm{\nabla }(p+\rho g \mathit{\text{z}})=\alpha v+\beta R_{\text{b}} (H_{\text{b}}-H),

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³), 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):

(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.

4.2.2 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 (long-term light rain and short-term heavy rain) were selected. Three typical soils (clay, silt, and sand) were analyzed in the test. Tab.3 shows the model parameters in the test, in which the rainfall intensity was disturbed up and down in steps of 20%. Tab.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%. Fig.12 shows the benchmark values of the soil moisture characteristic function used in the numerical tests.

4.2.3 Results and analyses

The numerical test considered the advancing distance of the wetting front as the model output. Fig.13 and Fig.14 show the test results. Within the range of parameter changes, the three parameters of rainfall intensity v, air-entry 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 ψ₀ 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). Fig.15 shows the calculation results. For the three typical soils under the two rainfall modes, the saturated hydraulic conductivity k_s, initial condition ψ₀, and air-entry pressure ψ_b are highly sensitive. The initial condition ψ₀ of clay is the most sensitive parameter under the long-term light rain mode. The initial condition ψ₀ 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 ψ₀ 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 ψ₀ 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 ψ₀, and air-entry value ψ_b should be ensured first in the numerical calculation.

4.3 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. Fig.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. Tab.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. Fig.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)). Fig.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 ψ₀, 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. Fig.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.

5 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 ψ₀ 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 ψ₀, 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.

References

[1] Park J S R, Cheung S W, Mai T. Multiscale simulations for multi-continuum Richards equations. Journal of Computational and Applied Mathematics, 2021, 397: 113648

[2] Richards L A. Capillary conduction of liquids through porous mediums. Physics, 1931, 1(5): 318–333

[3] Farthing M W, Ogden F L. Numerical solution of Richards’ equation: A review of advances and challenges. Soil Science Society of America Journal, 2017, 81(6): 1257–1269

[4] Li Z, Özgen-Xian I, Maina F Z. A mass-conservative predictor-corrector solution to the 1D Richards equation with adaptive time control. Journal of Hydrology (Amsterdam), 2021, 592: 125809

[5] Paniconi C, Putti M. Physically based modeling in catchment hydrology at 50: Survey and outlook. Water Resources Research, 2015, 51(9): 7090–7129

[6] Klute A. The determination of the hydraulic conductivity and diffusivity of unsaturated soils1. Soil Science, 1972, 113(4): 264–276

[7] Li X, Zhang L M, Fredlund D G. Wetting front advancing column test for measuring unsaturated hydraulic conductivity. Canadian Geotechnical Journal, 2009, 46(12): 1431–1445

[8] Liu Q, Xi P, Miao J, Li X, Wang K. Applicability of wetting front advancing method in the sand to silty clay soils. Soil and Foundation, 2020, 60(5): 1215–1225

[9] Ma D, Wang Q, Shao M. Analytical method for estimating soil hydraulic parameters from horizontal absorption. Soil Science Society of America Journal, 2009, 73(3): 727–736

[10] Su W, Cui Y J, Qin P J, Zhang F, Ye W M, Conil N. Application of instantaneous profile method to determine the hydraulic conductivity of unsaturated natural stiff clay. Engineering Geology, 2018, 243: 111–117

[11] Li X, Zhang Z, Zhang L, Zhang L, Wu L. Combining two methods for the measurement of hydraulic conductivity over a wide suction range. Computers and Geotechnics, 2021, 135: 104178

[12] Bruce R R, Klute A. The measurement of soil moisture diffusivity1. Soil Science Society of America Journal, 1956, 20(4): 458

[13] Morel-Seytoux H J, Khanji J. Prediction of imbibition in a horizontal column. Soil Science Society of America Journal, 1975, 39(4): 613–617

[14] NingLLikos W J. Unsaturated Soil Mechanics. New Jerse: Wiley, 2004

[15] BrooksR HCorey A T. Hydraulic Properties of Porous Media. Fort Collins: Hydrol Pap, 1964

[16] van Genuchten M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 1980, 44(5): 892–898

[17] Fredlund D G, Xing A. Equations for the soil-water characteristic curve. Canadian Geotechnical Journal, 1994, 31(4): 521–532

[18] Zha Y, Yang J, Shi L, Song X. Simulating one-dimensional unsaturated flow in heterogeneous soils with water content-based Richards equation. Vadose Zone Journal, 2013, 12(2): vzj2012.0109

[19] Tai B, Yue Z, Qi S, Wang P. Experimental and numerical investigation on thermal-moisture-mechanical behaviors on a new anti-frost cutting bed of high-speed railway in deep seasonally frozen ground regions under extreme climate. Computers and Geotechnics, 2021, 136: 104251

[20] Li X, Li X, Wu Y, Wu L, Yue Z. Selection criteria of mesh size and time step in fem analysis of highly nonlinear unsaturated seepage process. Computers and Geotechnics, 2022, 146: 104712

[21] Tai B, Wu Q, Zhang Z, Xu X. Study on thermal performance of novel asymmetric crushed-rock-based embankment on the qinghai-tibet railway in permafrost region. International Journal of Thermal Sciences, 2020, 152: 106333

[22] Fidkowski K J, Darmofal D L. Review of output-based error estimation and mesh adaptation in computational fluid dynamics. AIAA Journal, 2011, 49(4): 673–694

[23] Hénon P, Ramet P, Roman J. Pastix: A high-performance parallel direct solver for sparse symmetric positive definite systems. Parallel Computing, 2002, 28(2): 301–321

[24] LiX. An Overview of SuperLU. New York: ACM Transactions on Mathematical Software (TOMS), 2005

[25] List F, Radu F A. A study on iterative methods for solving Richards’ equation. Computational Geosciences, 2015, 20(2): 1–13

[26] Pettway J S, Schmidt J H, Stagg A K. Adaptive meshing in a mixed regime hydrologic simulation model. Computational Geosciences, 2010, 14(4): 665–674

[27] Pop I S, Radu F, Knabner P. Mixed finite elements for the Richards’ equation: Linearization procedure. Journal of Computational and Applied Mathematics, 2004, 168(1−2): 365–373

[28] Solin P, Kuraz M. Solving the nonstationary Richards equation with adaptive HP-FEM. Advances in Water Resources, 2011, 34(9): 1062–1081

[29] Valdez A R, Rocha B M, Chapiro G, Weber dos Santos R. Uncertainty quantification and sensitivity analysis for relative permeability models of two-phase flow in porous media. Journal of Petroleum Science Engineering, 2020, 192: 107297

[30] Sieber A, Uhlenbrook S. Sensitivity analyses of a distributed catchment model to verify the model structure. Journal of Hydrology (Amsterdam), 2005, 310(1−4): 216–235

[31] Berg S, Unsal E, Dijk H. Sensitivity and uncertainty analysis for parameterization of multiphase flow models. Transport in Porous Media, 2021, 140(1): 27–57

[32] Crosetto M, Tarantola S, Saltelli A. Sensitivity and uncertainty analysis in spatial modelling based on GIS. Agriculture, Ecosystems & Environment, 2000, 81(1): 71–79

[33] Morris M D. Factorial sampling plans for preliminary computational experiments. Technometrics, 1991, 33(2): 161–174

[34] Campolongo F, Tarantola S, Saltelli A. Tackling quantitatively large dimensionality problems. Computer Physics Communications, 1999, 117(1−2): 75–85

[35] Herman J D, Kollat J B, Reed P M, Wagener T. Technical note: Method of Morris effectively reduces the computational demands of global sensitivity analysis for distributed watershed models. Hydrology and Earth System Sciences, 2013, 17(7): 2893–2903

[36] Pang M, Xu R, Hu Z, Wang J, Wang Y. Uncertainty and sensitivity analysis of input conditions in a large shallow lake based on the latin hypercube sampling and morris methods. Water (Basel), 2021, 13(13): 1861

[37] King D M, Perera B J C. Morris method of sensitivity analysis applied to assess the importance of input variables on urban water supply yield—A case study. Journal of Hydrology (Amsterdam), 2013, 477: 17–32

[38] Ren J, Zhang W, Yang J. Morris sensitivity analysis for hydrothermal coupling parameters of embankment dam: A case study. Mathematical Problems in Engineering, 2019, 2019: 1–11

[39] Yi X, Zou R, Guo H. Global sensitivity analysis of a three-dimensional nutrients-algae dynamic model for a large shallow lake. Ecological Modelling, 2016, 327: 74–84

[40] Niessner J, Helmig R, Jakobs H, Roberts J E. Interface condition and linearization schemes in the newton iterations for two-phase flow in heterogeneous porous media. Advances in Water Resources, 2005, 28(7): 671–687

[41] Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 1976, 12(3): 513–522

[42] CrankJ. The Mathematics of Diffusion. New York: Oxford University Press, 1979

[43] Tracy F T. Clean two- and three-dimensional analytical solutions of Richards’ equation for testing numerical solvers. Water Resources Research, 2006, 42(8): W08503

[44] Gibson R E. A note on the constant head test to measure soil permeability in situ. Geotechnique, 1966, 16(3): 256–259

[45] Gibson R E. An extension to the theory of the constant head in situ permeability test. Geotechnique, 1970, 20(2): 193–197

[46] Powrie W. Soil Mechanics: Concepts and Applications. 2nd ed. Oxford: Taylor & Francis, 2013, 22(3): 381–382

[47] Saltelli A, Sobol I M. About the use of rank transformation in sensitivity analysis of model output. Reliability Engineering & System Safety, 1995, 50(3): 225–239

[48] Zádor J, Zsély I G, Turányi T. Local and global uncertainty analysis of complex chemical kinetic systems. Reliability Engineering & System Safety, 2006, 91(10-11): 1232–1240

[49] Francos A, Elorza F J, Bouraoui F, Bidoglio G, Galbiati L. Sensitivity analysis of distributed environmental simulation models: Understanding the model behaviour in hydrological studies at the catchment scale. Reliability Engineering & System Safety, 2003, 79(2): 205–218

[50] Chui T F M, Freyberg D L. Implementing hydrologic boundary conditions in a multiphysics model. Journal of Hydrologic Engineering, 2009, 14(12): 1374–1377

RIGHTS & PERMISSIONS

The Author(s) 2022. This article is published with open access at link.springer.com and journal.hep.com.cn