Finite analytic method for simulating water flow using water content-based Richards' equation

Abstract

Accurately simulating water flow movement in vadose zone is crucial for effective water resources assessment. Richards' equation, which describes the movement of water flow in the vadose zone, is highly nonlinear and challenging to solve. Existing numerical methods often face issues such as numerical dispersion, oscillation, and mass non-conservation when spatial and temporal discretization conditions are not appropriately configured. To address these problems and achieve accurate and stable numerical solutions, a finite analytic method based on water content-based Richards' equation (FAM-W) is proposed. The performance of the FAM-W is compared with analytical solutions, Finite Difference Method (FDM), and Finite Analytic Method based on the pressure Head-based Richards' equation (FAM-H). Compared to analytical solution and other numerical methods (FDM and FAM-H), FAM-W demonstrates superior accuracy and efficiency in controlling mass balance errors, regardless of spatial step sizes. This study introduces a novel approach for modelling water flow in the vadose zone, offering significant benefits for water resources management.

Graphical abstract

Keywords

Finite analytic method / Vadose zone / Soil moisture / Finite difference method / Analytical solution / Richards' equation / Water resources management

Introduction

The vadose zone is a critical component of the Earth's hydrological cycle (Vereecken et al. 2022). Extending from land surface to water table, this zone plays a pivotal role in various environmental processes, including groundwater recharge, evapotranspiration, pollutant transport, and ecosystem functioning (Li et al. 2019). Accurately simulation of water flow is essential for effective water resource management and environmental protection (Tu et al. 2020; Wang et al. 2021). However, solving Richards' equation, which describes water movement in the vadose zone, presents significant challenges due to its highly nonlinear nature (Berardi et al. 2023b; Sun et al. 2024; Zhang et al. 2020; Zhang et al. 2018). The nonlinear parameters of Richards' equation further complicate its solution (Berardi et al. 2023a), making the analytical methods impractical and limiting their applicability to specific, simplified scenarios (Timsina, 2024). As a result, advanced numerical methods are required to address these nonlinearities in both the equation and associated parameter estimation (Younes et al. 2022; Ren et al. 2021; Wang et al. 2019).

Richards' equation can be expressed in three different forms: Pressure head-form, mixed-form, and water content-form (Berardi and Difonzo, 2022). Each form has advantages and disadvantages depending on simulation scenarios. The widely used pressure head-form is suitable for simulating water flow in soils with varying saturation levels, from fully saturated to partially saturated. However, Zha et al. (2019) noted that this form significantly reduces computational efficiency when applied to dry soils with low water content. The pressure-head, which represents the negative pressure exerted by the soil matrix on soil water due to capillary and adsorptive forces (Haverkamp et al. 2016), often leads to significant mass balance errors due to the highly nonlinear of the soil water retention curve (Zha et al. 2017). The mixed-form of Richards' equation better controls mass balance errors, as demonstrated by Zhang et al. (2016). Suk and Park (2019) enhanced the mixed-form by employing a Picard iteration method alongside finite difference and finite element methods, effectively mitigating the mass balance errors. Nevertheless, Zha et al. (2013) identified significant computational errors when using the mixed-form for infiltration processes in dry soils. The water content-based Richards' equation is frequently used in meteorology due to its effectiveness in controlling mass balance errors, even with larger spatial step sizes (Zeng and Decker, 2009). It also exhibits high numerical stability (Shen and Phanikumar, 2010) and reduces computational time, making it advantageous for large-scale meteorological modelling.

Currently, the Finite Element Method (FEM) and Finite Difference Method (FDM) are the most commonly used numerical methods for solving Richards' equation (Zhang et al. 2015). FEM employs shape functions to discretize the domain, offering flexibility in handling complex geometries. FDM, which directly discretizes the equation using a grid, is straightforward and easy to implement for regular grids. However, both methods have inherent limitations, including numerical dispersion, numerical oscillation, and mass non-conservation when inappropriate spatial and temporal discretization conditions are applied (Celia et al. 1990). To address these issues, researchers have developed improved numerical methods (Kirkland et al, 1992; Suk and Park, 2019; Szymkiewicz, 2009). In the 1980s, Chen and Chen (1984) introduced the finite analytic method (FAM), which preserves the physical characteristics of the original problem and provides monotonic, oscillation-free, and robust stability (Civan, 2009). Zhang et al. (2015) applied FAM to the pressure head form (FAM-H), demonstrating superior accuracy compared to FDM. However, as a non-mass conservation method, FAM-H does not guarantee mass balance (Tsai et al. 1993). Zhang et al. (2016) extended FAM to the mixed-form (FAM-M) and yielded the highest accuracy among FDM and FAM-H. Nevertheless, FAM-M introduce errors when simulating soil moisture in dry conditions (Zhang et al. 2021).

Given the advantages of the water content-based Richards' equation in representing infiltration processes, no study has applied FAM to this form. This study aims to fill that gap by developing a Finite Analytic Method (FAM) for the water content-based Richards' equation (FAM-W). The study first derives the finite analytical computational format for the equation and evaluates its accuracy against analytical solutions. Since FDM and FEM yield comparable accuracy in simulating soil water content movement in the vadose zone (Gottardi and Venutelli, 1993), FEM is not employed here. Instead, numerical solutions from FAM-W are compared with those from FDM and FAM-H. Local and mass balance errors of the three numerical methods are also quantitatively analysed. This research introduces a novel approach for simulating water movement in the vadose zone, offering significant benefits for water resource management.

1 Description of the model

1.1 Derivation of FAM-W

The governing equation for the water content-based Richards' equation is expressed as:

(1)$ \frac{{\partial \theta }}{{\partial t}} = \frac{\partial }{{\partial {\textit{z}}}}\left[ {D\left( \theta \right)\frac{{\partial \theta }}{{\partial {\textit{z}}}}} \right] + \frac{{\partial K\left( \theta \right)}}{{\partial {\textit{z}}}} $

Where: $ \theta $ represents the soil water content (cm³/cm³), $ D\left( \theta \right) = K\left( \theta \right)\dfrac{{\partial h}}{{\partial \theta }} $ is the hydraulic diffusivity (cm²/h), and $ K\left( \theta \right) $ is the unsaturated hydraulic conductivity (cm/h), which used the Gardner model to describe (Gardner, 1958), $ {\textit{z}} $ denotes the vertical dimension with the upward direction designated as positive. Further, the unsaturated hydraulic conductivity $ K\left( \theta \right) $ is given by:

(2)$ K\left( \theta \right) = {K_s}\frac{{(\theta - {\theta _r})}}{{({\theta _s} - {\theta _r})}} $

(3)$ \theta = exp\left( {\alpha h} \right)\left( {{\theta _s} - {\theta _r}} \right) + {\theta _r} $

Where: $ {K_{{s}}} $ is the saturated hydraulic conductivity (cm/h), and $ {\theta _s} $ and $ {\theta _r} $ are the saturated and residual soil water content (cm³/cm³), respectively. $ h $ is the pressure head (cm), and $ \alpha $ is a soil index parameter related to the pore-size distribution (1/cm).

(4)$ D\left( \theta \right) = \frac{{K\left( h \right)}}{{{{d\theta } \mathord{\left/ {\vphantom {{d\theta } {dh}}} \right. } {dh}}}} = \frac{{{K_s}}}{{\alpha \left( {{\theta _s} - {\theta _r}} \right)}} $

(5)$ \frac{{\partial K}}{{\partial {\textit{z}}}} = \frac{{{K_s}}}{{\left( {{\theta _s} - {\theta _r}} \right)}}\left( {\frac{{\partial \theta }}{{\partial {\textit{z}}}}} \right) $

As a consequence, Eq. (1) can be expressed as:

(6)$ \frac{{{\partial ^2}\theta }}{{\partial {{\textit{z}}^2}}} + {x_1}\frac{{\partial \theta }}{{\partial {\textit{z}}}} = {x_2}\frac{{\partial \theta }}{{\partial t}} $

Where: $ {x_1} = \alpha $ and $ \;{x_2} = \dfrac{{{\theta _s} - {\theta _r}}}{{{K_s}}} $. Eq. (6) can be expressed as:

(7)$ {x_2}\frac{{(\theta _i^n - \theta _i^{n - 1})}}{{\Delta t}} \approx {\left(\frac{{{\partial ^2}\theta }}{{\partial {{\textit{z}}^2}}}\right)^n} + {x_1}{\left(\frac{{\partial \theta }}{{\partial {\textit{z}}}}\right)^n} $

Where:

(8)$ f = {x_2}\frac{{(\theta _i^n - \theta _i^{n - 1})}}{{\Delta t}} $

(9)$ B = - \frac{{{x_1}}}{2} $

Eq. (7) can be expressed as:

(10)$ {\left( {\frac{{{\partial ^2}\theta }}{{\partial {{\textit{z}}^2}}}} \right)^n} - 2B{\left( {\frac{{\partial \theta }}{{\partial {\textit{z}}}}} \right)^n} = f $

By discretizing the equation for a local element as illustrated in Fig. 1 (modified from Zhang et al. 2015), the equation becomes:

(11)$ {\theta ^n}\left( {\textit{z}} \right) = {b_1} + {b_2}{e^{2B{\textit{z}}}} - \frac{f}{{2B}}{\textit{z}} $

Constants $ {b_1} $ and $ {b_2} $ in Eq. (12) can be determined for each local element as:

(12)$ \begin{split} & {b}_{1}=\frac{{\theta }^{n}\left(i+1\right)+{\theta }^{n}(i-1)}{2}- \\& ch(2B\Delta {\textit{z}})\left[\frac{{\theta }^{n}(i+1)-{\theta }^{n}(i-1)}{2sh(2B\Delta {\textit{z}})}+\frac{f\Delta {\textit{z}}}{2sh(2B\Delta {\textit{z}})}\right]\end{split} $

(13)$ {b_2} = \frac{{{\theta ^n}(i + 1) - {\theta ^n}(i - 1)}}{{2sh(2B\Delta {\textit{z}})}} + \frac{{f\Delta {\textit{z}}}}{{2sh(2B\Delta {\textit{z}})}} $

According to Zhang et al. (2015), the following computational format can be obtained:

(14)$ \begin{split} {\theta ^n}(i) = & \frac{{[1 - th(B\Delta {\textit{z}})]}}{{\gamma + 2}}{\theta ^n}(i + 1) +\\& \frac{{[1 + th(B\Delta {\textit{z}})]}}{{\gamma + 2}}{\theta ^n}(i - 1) + \frac{\gamma }{{\gamma + 2}}{\theta ^{n - 1}}(i) \end{split}$

Where:

(15)$ \gamma = \frac{{\Delta {\textit{z}}th(B\Delta {\textit{z}})}}{{B\Delta t}}{x_2} $

1.2 FAM-H and FDM

To evaluate the performance of FAM-W, we compare it with the results from FAM-H (Zhang et al. 2015) and FDM (Celia, 1990). First, we derive the computational format of FAM-H, The equation for FAM-H is expressed as:

(16)$ C\left( h \right)\frac{{\partial h}}{{\partial t}} = \frac{\partial }{{\partial {\textit{z}}}}\left[ {K\left( h \right)\frac{{\partial h}}{{\partial {\textit{z}}}}} \right] + \frac{{\partial K\left( h \right)}}{{\partial {\textit{z}}}} $

Where: $ C\left( h \right) = d\theta /dh $ represents the specific moisture capacity function (1/cm), $h$ is the pressure head (cm), $ K(h) $ is the unsaturated hydraulic conductivity (cm/h).

(17)$ K\left( h \right) = {K_s}\exp (\alpha h) $

Using Kirchhoff transformation, Eq. (16) becomes:

(18)$ \frac{{C(h)}}{{K(h)}}\frac{{\partial u}}{{\partial t}} = \frac{{{\partial ^2}u}}{{\partial {{\textit{z}}^2}}} + \frac{{dK(h)}}{{dh}}\frac{1}{{K(h)}}\frac{{\partial u}}{{\partial {\textit{z}}}} $

According to Zhang et al. (2015), the following computational format for FAM-H is obtained as:

(19)$ \begin{split} {u^n}(i) =& \frac{{[1 - th(B\Delta {\textit{z}})]}}{{2 + \lambda }}{u^n}(i + 1) + \\& \frac{{[1 + th(B\Delta {\textit{z}})]}}{{2 + \lambda }}{u^n}(i - 1) + \frac{\lambda }{{2 + \lambda }}{u^{n - 1}}(i) \end{split}$

Where:

(20)$ \lambda = \frac{{\Delta {\textit{z}}th(B\Delta {\textit{z}})}}{{B\Delta t}}\frac{{{C^{n - 1}}(h)}}{{{K^{n - 1}}(h)}} $

(21)$ B = - \frac{1}{{2K(h)}}\frac{{dK(h)}}{{dh}} $

Boundary conditions are incorporated using a forward difference scheme with second-order approximation (Eq. 22).

(22)$ \frac{{\partial \theta }}{{\partial {\textit{z}}}} = \frac{1}{{2\Delta {\textit{z}}}}[\theta (i - 2) - 4\theta (i - 1) + 3\theta (i)] $

Each node within the computational domain obtained the computational format of Eqs. (14) and (19), along with the boundary conditions, to form a system of equations. The numerical solutions $ {\theta }^{n}\left(i\right) $ can be obtained using the Successive Over-Relaxation (SOR) method, which iterates until convergence. For a detailed derivation of FDM, please refer to the paper by Celia et al. (1990).

1.3 Evaluation of the model

To evaluate the computational accuracy of FAM-W, it was compared to the results of analytical solutions which serve as reference values. A number of studies have derived analytical solutions for simulating soil moisture in the vadose zone. Among these, Srivastava and Yeh (1991) proposed a widely-used solution for vertical infiltration in homogeneous, one-dimensional soil. This study employs their analytical solution as the benchmark for evaluating these numerical methods. For the analytical solution, the soil is homogeneous and its properties are defined as follows: Thickness = 1 m, bulk density = 1.4 g/cm³, porosity = 0.4, saturated hydraulic conductivity = 1 cm/h, and $\alpha $= 0.06 cm⁻¹. The flux at the upper boundary is initially set to 0.1 cm/h. Once the soil profile reaches a steady state, we take the soil water content along with the soil profile as the initial condition. Then, the flux at the upper boundary is changed to be 0.9 cm/h, with the lower boundary saturated. The numerical performances of FAM-H, FAM-M, and FDM are evaluated under different spatial step conditions ($\Delta {\textit{z}}$=1, 2, and 5 cm), with results compared to with the analytical solution. Local errors (Eq. 23) and mass balance errors (Eq. 24) are calculated to analyze the accuracy of each method.

(23)$ \delta \left( {{\textit{z}},t} \right) = \left( {1 - \frac{{\theta {{({\textit{z}},t)}_{numerical}}}}{{\theta {{({\textit{z}},t)}_{analytical}}}}} \right) \times 100\% $

(24)$ \varepsilon \left( t \right) = \left[1 - \frac{{\left( {M_{numerical}^t - M_{numerical}^0} \right)}}{{\left( {M_{analytical}^t - M_{analytical}^0} \right)}}\right] \times 100\% $

Where: $ \theta {({\textit{z}},t)_{numerical}} $ and $ \theta {({\textit{z}},t)_{analytical}} $ represent the soil water content from the numerical model and analytical solution, respectively. $ M_{numerical}^0 $ and $ M_{numerical}^t $ are the initial mass and the mass at the time $t$ of numerical results in the domain, respectively. $ M_{analytical}^0 $ and $ M_{analytical}^t $ are the initial mass and the mass at the time $t$ of analytical solution in the domain, respectively.

2 Results

2.1 Comparison with analytical solution and the results of FAM-M, FAM-H, and FDM

Fig. 2 displays the numerical solutions obtained from FAM-W, FAM-H, FDM, and the analytical solution. As shown in Fig. 2(a), when the spatial step size is set to 1 cm, all three methods yield satisfactory numerical results. However, when the spatial step size increased to 2 cm, the results of FAM-W are more accurate than those of FAM-H and FDM. This trend becomes more pronounced as the spatial step size is further increased to 5 cm. These findings indicate that FAM-W is less sensitive to the changes in spatial step size compared to FAM-H and FDM. At a spatial step size of 2 cm, the results of FAM-H tend to overestimate the analytical solution, while FDM underestimates it (Fig. 2(b)). At 5 cm, the overestimation by FAM-H and underestimation by FDM are more evident (Fig. 2(c)). Moreover, the errors of FDM are significantly larger than those of FAM-H, indicating that FDM is more sensitive to variations in the spatial step size.

2.2 The local errors of FAM-M, FAM-H, and FDM

Figs. 3, 4, and 5 provide a quantitative comparison of local errors in the soil profiles obtained from FAM-W, FAM-H, and FDM.

FAM-W: As seen in Fig. 3(a), the maximum local error of FAM-W does not exceed 1.5%. The errors primarily occur during the initial stages of infiltration but decrease rapidly over time, indicating that FAM-W produces highly accurate numerical solutions. All local errors are positive, suggesting that FAM-W slightly underestimates the analytical solutions. The maximum local error occurs at 91 cm at t=1 h, at 63 cm at t=5 h, at 21 cm at t=20 h, and at 1 cm at t=100 h. This shows that the largest local error is typically found at the wetting front of the infiltration. Figs. 3(b) and 3(c) show similar distributions of local error. As the spatial step size increases, the maximum local errors also increase. For example, at the spatial step size of 5 cm, the maximum local error reaches 6.1%.

FAM-H: Fig. 4(a) shows that FAM-H results in considerable local errors during the initial stages of infiltration process, which decrease over time. The positive local errors suggest that FAM-H underestimates the soil moisture compared to those of the analytical solution. As illustrated in Fig. 4(a), the maximum local error did not exceed 3.5% during the initial stage. As shown in Figs. 4(b) and 4(c), the maximum local errors increase as the spatial step size grows. For example, at the spatial step size of 5 cm, the maximum local error reaches 20.0%.

FDM: Fig. 5(a) demonstrates significant local errors in FDM at the initial stages of infiltration, with these errors gradually decreasing over time. The negative local errors indicate that FDM overestimates the analytical solution. The largest local error shown in Fig. 5(a) exceeds 6% that occurs in the initial infiltration stage. As shown in Figs. 5(b) and 5(c), local errors increase with larger spatial step sizes. At 5 cm, the absolute value of the maximum local error exceeds 30%. These results highlight that FDM performs satisfactorily only when smaller spatial step sizes are used. The above-mentioned results indicate that, compared to FAM-H and FDM, FAM-W demonstrates superior control over local errors and yields more accurate numerical solutions.

2.3 The mass balance errors of FAM-W and FAM-H

To further evaluate the performance of FAM-W, the mass balance errors are also calculated. Since FDM has the largest local errors compared to FAM-M and FAM-H, we focus on comparing the mass balance errors between FAM-W and FAM-H. The results of mass balance errors for FAM-W and FAM-H are shown in Fig. 6. For FAM-W, Fig. 6a shows the mass balance errors decrease over time. The maximum absolute mass balance error does not exceed 1.5% when the spatial step size is 1 cm. As the spatial step size increases, the maximum mass balance error increases, reaching 6% at a spatial step size of 5 cm. The mass balance errors for FAM-H follow a similar trend but are generally larger than those of FAM-W. The maximum mass balance errors for FAM-W are 1.23%, 1.88%, and 5.11% at spatial step sizes of 1 cm, 2 cm, and 5 cm, respectively. In contrast, the maximum mass balance errors for FAM-H are 1.90%, 4.67%, and 13.79% for grid sizes of 1 cm, 2 cm, and 5 cm, respectively. These results suggest that FAM-W has better control over mass balance errors than FAM-H.

3 Discussion

Examining the local errors, it is observed that the maximum local error across the three methods occurs at the wetting front. This finding indicates that all methods experience errors in regions where there are significant gradient variations during the infiltration process. According to Kumar et al. (2023), the accuracy of FDM in simulating soil moisture dynamics is constrained by its nonlinear character. FDM solves differential equations by approximating derivatives with finite differences. However, the error associated with employing finite differences to calculate $\dfrac{{\partial K\left( h \right)}}{{\partial {\textit{z}}}}$ proves to be a limitation of FDM. The accuracy of FDM is highly sensitive to grid size, with coarse grids resulting in greater discrepancies. When the grid is too coarse, truncation errors increase because the derivative approximation becomes less precise. For instance, FDM struggles to maintain low local errors, as shown in Fig. 5, where the maximum local error exceeds 30% at a spatial step size of 5 cm. Therefore, FDM can only yield satisfactory results when both time and spatial steps are sufficiently small, as shown in Fig. 2a.

FAM-H, which uses the Kirchhoff transform, does not involve the computation of the $\dfrac{{\partial K\left( h \right)}}{{\partial {\textit{z}}}}$, making it less sensitive to spatial step size than FDM. However, FAM-H generally has larger local and mass balance errors compared to FAM-W. This is primarily due to the high non-linearity of $C(h)$ and $K\left( h \right)$ in the Richards' equation (Zha et al. 2013) and the use of finite difference method for calculating the $C(h)\dfrac{{\partial h}}{{\partial t}}$ term in FAM-H, which leads to large truncation errors when larger spatial steps are used (Zhang et al. 2016). On the other hand, FAM-W provides the highest accuracy in numerical solutions and better control over mass balance errors, especially for larger grid sizes. This is due to the fact that the water content-based Richards' equation is inherently mass-conservative and less non-linear in terms of inter-nodal averaged hydraulic diffusivity, which reduces numerical difficulties (Zeng et al. 2018).

4 Conclusion

In this study, we developed a finite analytic method based on the water content-based Richards' equation (FAM-W) to simulate water flow in the vadose zone. This new approach enriches the theory of FAM. A comparative analysis was conducted to evaluate the performance of FAM-W in relation to analytical solutions and numerical solutions (FAM-H and FDM). All three methods exhibit errors when simulating the infiltration process at the wetting fronts. However, FAM-W delivers the most accurate numerical solutions for simulating the movement of the water flow in the vadose zone. The performance of FAM-H is the second-best, while FDM shows the poorest performance. FAM-W demonstrates the least sensitivity to temporal and spatial step sizes, followed by FAM-H, while the FDM method is the most sensitive to the step size. In addition, the FAM-W method offers better control over mass balance error.

References

[1] Berardi M, Difonzo FV. 2022. A quadrature-based scheme for numerical solutions to Kirchhoff transformed Richards' equation. Journal of Computational Dynamics, 9: 69−84. DOI: 10.3934/jcd.2022001.

[2] Berardi M, Difonzo FV, Guglielmi R. 2023a. A preliminary model for optimal control of moisture content in unsaturated soils. Computational Geosciences, 27: 1133−1144. DOI: 10.1007/s10596-023-10250-1.

[3] Berardi M, Difonzo FV, Pellegrino SF. 2023b. A numerical method for a nonlocal form of Richards' equation based on peridynamic theory. Computers & Mathematics with Applications, 143: 23−32. DOI: 10.1016/j.camwa.2023.04.032.

[4] Celia MA, Bouloutas ET, Zarba RL. 1990. A general mass‐conservative numerical solution for the unsaturated flow equation. Water resources research, 26: 1483−1496. DOI: 10.1029/WR026i007p01483.

[5] Chen CJ, Chen HC. 1984. Finite analytic numerical method for unsteady two-dimensional Navier-Stokes equations. Journal of Computational Physics, 53: 209−226. DOI: 10.1016/0021-9991(84)90038-X.

[6] Civan F. 2009. Practical finite-analytic method for solving differential equations by compact numerical schemes. Numerical Methods for Partial Differential Equations, 25: 347−379. DOI: 10.1002/num.20346.

[7] Gardner W. 1958. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Science, 85: 228−232. DOI: 10.1097/00010694-195804000-00006.

[8] Gottardi G, Venutelli M. 1993. Richards: Computer program for the numerical simulation of one-dimensional infiltration into unsaturated soil. Computers & Geosciences, 19: 1239−1266. DOI: 10.1016/0098-3004(93)90028-4.

[9] Haverkamp R, Debionne S, Angulo-Jaramillo R, et al. 2016. Soil properties and moisture movement in the unsaturated zone. In "The handbook of groundwater engineering", 167-208. CRC Press. DOI: 10.1201/9781420006001.

[10] Kirkland MR, Hills R, Wierenga P. 1992. Algorithms for solving Richards' equation for variably saturated soils. Water Resources Research, 28: 2049−2058. DOI: 10.1029/92WR00802.

[11] Kumar R, Tiwari R, Prasad R. 2023. Numerical solution of partial differential equations: Finite difference method. Mathematics and Computer Science, 1: 353−372. DOI: 10.1007/978-1-4899-7278-1.

[12] Li M, Li J, Singh VP, et al. 2019. Efficient allocation of agricultural land and water resources for soil environment protection using a mixed optimization-simulation approach under uncertainty. Geoderma, 353: 55−69. DOI: 10.1016/j.geoderma.2019.06.023.

[13] Ren S, ZhangY, Jim Yeh TC, et al. 2021. Multiscale hydraulic conductivity characterization in a fractured granitic aquifer: The evaluation of scale effect. Water Resources Research, 57: e2020WR028482. DOI: 10.1029/2020WR028482.

[14] Shen C, Phanikumar MS. 2010. A process-based, distributed hydrologic model based on a large-scale method for surface–subsurface coupling. Advances in Water Resources, 33: 1524−1541. DOI: 10.1016/j.advwatres.2010.09.002.

[15] Srivastava R, Yeh TCJ. 1991. Analytical solutions for one‐dimensional, transient infiltration toward the water table in homogeneous and layered soils. Water Resources Research, 27: 753−762. DOI: 10.1029/90WR02772.

[16] Suk H, Park E. 2019. Numerical solution of the Kirchhoff-transformed Richards equation for simulating variably saturated flow in heterogeneous layered porous media. Journal of Hydrology, 579: 124213. DOI: 10.1016/j.jhydrol.2019.124213.

[17] Sun Y, Li D, Jiao L, et al. 2024. Simulation of saturated–unsaturated seepage problems via the virtual element method. Computers and Geotechnics, 171: 106326. DOI: 10.1016/j.compgeo.2024.106326.

[18] Szymkiewicz A. 2009. Approximation of internodal conductivities in numerical simulation of one‐dimensional infiltration, drainage, and capillary rise in unsaturated soils. Water Resources Research, 45(10): W10403. DOI: 10.1029/2008WR007654.

[19] Timsina RC. 2024. Infiltration-Induced landslide: An application of richards equation. Journal of Nepal Mathematical Society, 7: 86−99. DOI: 10.3126/jnms.v7i1.67490.

[20] Tsai WF, Chen CJ, Tien HC. 1993. Finite analytic numerical solutions for unsaturated flow with irregular boundaries. Journal of Hydraulic Engineering, 119: 1274−1298. DOI: 10.1061/(ASCE)0733-9429(1993)119:11(1274).

[21] Tu MC, Wadzuk B, Traver R. 2020. Methodology to simulate unsaturated zone hydrology in Storm Water Management Model (SWMM) for green infrastructure design and evaluation. PLoS One, 15: e0235528. DOI: 10.1371/journal.pone.0235528.

[22] Vereecken H, Amelung W, Bauke SL, et al. 2022. Soil hydrology in the Earth system. Nature Reviews Earth & Environment, 3: 573−587. DOI: 10.1038/s43017-022-00324-6.

[23] Wang YL, Yeh TCJ, Wen JC, et al. 2019. Resolution and ergodicity issues of river stage tomography with different excitations. Water Resources Research, 55: 4974−4993. DOI: 10.1029/2018WR023204.

[24] Wang YL, Yeh TCJ, Xu D, et al. 2021. Stochastic analysis of oscillatory hydraulic tomography. Journal of Hydrology, 596: 126105. DOI: 10.1016/j.jhydrol.2021.126105.

[25] Younes A, Koohbor B, Belfort B, et al. 2022. Modeling variable-density flow in saturated-unsaturated porous media: An advanced numerical model. Advances in Water Resources, 159: 104077. DOI: 10.1016/j.advwatres.2021.104077.

[26] Zeng J, Zha Y, Yang J. 2018. Switching the Richards' equation for modeling soil water movement under unfavorable conditions. Journal of Hydrology, 563: 942−949. DOI: 10.1016/j.jhydrol.2018.06.069.

[27] Zeng X, Decker M. 2009. Improving the numerical solution of soil moisture–based Richards equation for land models with a deep or shallow water table. Journal of Hydrometeorology, 10: 308−319. DOI: 10.1175/2008JHM1011.1.

[28] Zha Y, Shi L, Ye M, et al. 2013. A generalized Ross method for two-and three-dimensional variably saturated flow. Advances in water resources, 54: 67−77. DOI: 10.1016/j.advwatres.2013.01.002.

[29] Zha Y, Yang J, Yin L, et al. 2017. A modified Picard iteration scheme for overcoming numerical difficulties of simulating infiltration into dry soil. Journal of hydrology, 551: 56−69. DOI: 10.1016/j.jhydrol.2017.05.053.

[30] Zha Y, Yang J, Zeng J, et al. 2019. Review of numerical solution of Richardson–Richards equation for variably saturated flow in soils. Wiley Interdisciplinary Reviews: Water, 6: e1364. DOI: 10.1002/wat2.1364.

[31] Zhang Z, Wang W, Chen L, et al. 2015. Finite analytic method for solving the unsaturated flow equation. Vadose Zone Journal, 14: 1-10. DOI: 10.2136/vzj2014.06.0073.

[32] Zhang Z, Wang W, Gong C, et al. 2020. Finite analytic method: Analysis of one-dimensional vertical unsaturated flow in layered soils. Journal of Hydrology, 597: 125716. DOI: 10.1016/j.jhydrol.2020.125716.

[33] Zhang Z, Wang W, Gong C, et al. 2018. Finite analytic method for modeling variably saturated flows. Science of the Total Environment, 621: 1151−1162. DOI: 10.1016/j.scitotenv.2017.10.112.

[34] Zhang Z, Wang W, Lu Y. 2021. Improved finite analytic method to simulate soil water movement in vadose zones. Transactions of the Chinese Society of Agricultural Engineering, 37: 55−61. DOI: 10.11975/j.issn.1002-6819.2021.18.007.

[35] Zhang Z, Wang W, Yeh TCJ, et al. 2016. Finite analytic method based on mixed-form Richards' equation for simulating water flow in vadose zone. Journal of Hydrology, 537: 146−156. DOI: 10.1016/j.jhydrol.2016.03.035.

RIGHTS & PERMISSIONS

Journal of Groundwater Science and Engineering Editorial Office