Nonlinear Richards Equation Research Articles

Localized RBF methods for modeling infiltration using the Kirchhoff-transformed Richards equation

We develop a new approach to solve the nonlinear Richards equation based on the Kirchhoff transformation and localized radial basis function (LRBF) techniques. Our aim is to reduce the nonlinearity of the governing equation and apply LRBF methods for modeling unsaturated flow through heterogeneous soils. In our methodology, we propose special techniques which deal with the heterogeneity of the medium in order to apply the Kirchhoff transformation where we used the Brooks and Corey model for the capillary pressure function and a power-law relation in saturation for the relative permeability function. The new approach allows us to avoid the technical issues encountered in the Kirchhoff transformation due to soil heterogeneity in order to reduce the nonlinearity of the model equation. The resulting Kirchhoff-transformed Richards equation is solved using LRBF methods which have advantages in terms of computational cost since they do not require mesh generation. Furthermore, LRBF techniques lead to a system with a sparse matrix which allows us to avoid ill-conditioned issues. To validate the developed approach for predicting the dynamics of unsaturated flow in porous media, numerical experiments are performed in one, two, and three-dimensional soils. The numerical results demonstrate the efficiency and accuracy of the proposed techniques for modeling infiltration through heterogeneous soils.

Prediction of numerical homogenization using deep learning for the Richards equation

For the nonlinear Richards equation as an unsaturated flow through heterogeneous media, we build a new coarse-scale approximation algorithm utilizing numerical homogenization. This approach follows deep neural networks (DNNs) to quickly and frequently calculate macroscopic parameters. More specifically, we train neural networks with a training set consisting of stochastic permeability realizations and corresponding computed macroscopic targets (effective permeability tensor, homogenized stiffness matrix, and right-hand side vector). Our proposed deep learning scheme develops nonlinear maps between such permeability fields and macroscopic characteristics, and the treatment for Richards equation’s nonlinearity is included in the predicted coarse-scale homogenized stiffness matrix, which is a novelty. This strategy’s good performance is demonstrated by several numerical tests in two-dimensional model problems, for predictions of the macroscopic properties and consequently solutions.

A Novel Radial Basis Function Approach for Infiltration-Induced Landslides in Unsaturated Soils

In this article, the modeling of infiltration--induced landslides, in unsaturated soils using the radial basis function (RBF) method, is presented. A novel approach based on the RBF method is proposed to deal with the nonlinear hydrological process in the unsaturated zone. The RBF is first adopted for curve fitting to build the representation of the soil water characteristic curve (SWCC) that corresponds to the best estimate of the relationship between volumetric water content and matric suction. The meshless method with the RBF is then applied to solve the nonlinear Richards equation with the infiltration boundary conditions. Additionally, the fictitious time integration method is adopted in the meshless method with the RBF for tackling the nonlinearity. To model the stability of the landslide, the stability analysis of infinite slope coupled with the nonlinear Richards equation considering the fluctuation of transient pore water pressure is developed. The validation of the proposed approach is accomplished by comparing with exact solutions. The comparative analysis of the factor of safety using the Gardner model, the van Genuchten model and the proposed RBF model is provided. Results illustrate that the RBF is advantageous for reconstructing the SWCC with better estimation of the relationship than conventional parametric Gardner and van Genuchten models. We also found that the computed safety factors significantly depend on the representation of the SWCC. Finally, the stability of landslides is highly affected by matric potential in unsaturated soils during the infiltration process.

МЕТОД IДЕНТИФIКАЦIЇ ПОТУЖНОСТI ДЖЕРЕЛ В РIВНЯННI РIЧАРДСА

In this paper a one-dimensional nonlinear Richards equation describing fluid flow in porous medium with inserted equalpowered sources is studied. An experimental iterational method is proposed to find source power to minimize the deviation of received humidity values from target values. Modeling was performed using numerical difference approximation of derivatives, resulting into a system of nonlinear equations with dependence from previous time step. The offered method allows to perform modeling for different source power values, and chooses the most suitable one.Iterations stop when they reach average modular difference value less than calculation error of numerical difference scheme. Here explicit scheme was used to save time, equations were tested for unsaturated medium only to avoid flooding the area, so source power is tested with given limitations. Results of simulations and choice for next source power approximations are described and compared until solution is found. This approach is considered as experimental so we plan to perform more analysis in the future.

A Solution of Richards’ Equation by Generalized Finite Differences for Stationary Flow in a Dam

The accurate description of the flow of water in porous media is of the greatest importance due to its numerous applications in several areas (groundwater, soil mechanics, etc.). The nonlinear Richards equation is often used as the governing equation that describes this phenomenon and a large number of research studies aimed to solve it numerically. However, due to the nonlinearity of the constitutive expressions for permeability, it remains a challenging modeling problem. In this paper, the stationary form of Richards’ equation used in saturated soils is solved by two numerical methods: generalized finite differences, an emerging method that has been successfully applied to the transient case, and a finite element method, for benchmarking. The nonlinearity of the solution in both cases is handled using a Newtonian iteration. The comparative results show that a generalized finite difference iteration yields satisfactory results in a standard test problem with a singularity at the boundary.

Algebraic model for one-dimensional horizontal water flow with arbitrary initial soil water content

A simple analytical solution of the equation for the one-dimensional horizontal permeability of soil water is important for estimating the hydraulic properties of soil. Our main objective was to develop analytical solutions to the nonlinear Richards equation, with constant-saturation upper boundary and an arbitrary initial soil water content (SWC) for horizontal absorption. We estimated the infiltration rate based on the hypothesis proposed by Parlange and carried out a series of transformations based on the Brooks–Corey model to obtain a theoretical function of the one-dimensional movement of water in unsaturated soil under an arbitrary initial SWC. The algebraic analytical solutions were simple, and the selection of the initial SWC was arbitrary. We assumed three scenarios of linear distributions of initial SWC, and Hydrus-1D software was used to simulate horizontal infiltration. Based on the comparison of algebraic and numerical results, the corrected algebraic model was proposed and verified by the arbitrary initial water content conditions when the maximum SWC was less than the half of saturated water content. The proposed method provides a description of horizontal infiltration for the heterogeneous initial SWCs.

On transport through heterogeneous media: application of conjugated reciprocal transformations

Conjugation of reciprocal transformations is used to solve a class of boundary value problems involving a source term relevant to water transport through a heterogeneous medium with a volumetric extraction mechanism. The main emphasis is the solution method that involves conjugation of reciprocal transformations, as well as other changes of variable, applied to a newly identified integrable model. The transport equation is a version of the nonlinear Richards equation, based on Buckingham’s extension of Darcy’s law to unsaturated media, now allowing for heterogeneous transport coefficients. In addition, there is a plant-root water extraction term that depends on both water content and position. The nonlinear boundary conditions have prescribed flux at one boundary and zero flux at the other barrier boundary of a one-dimensional heterogeneous medium.

A control volume scheme using compact integrated radial basis function stencils for solving the Richards equation

A new control volume approach is developed based on compact integrated radial basis function (CIRBF) stencils for solution of the highly nonlinear Richards equation describing transient water flow in variably saturated soils. Unlike the conventional control volume method, which is regarded as second-order accurate, the proposed approach has high-order accuracy owing to the use of a compact integrated radial basis function approximation that enables improved flux predictions. The method is used to solve the Richards equation for transient flow in 1D homogeneous and heterogeneous soil profiles. Numerical results for different boundary conditions, initial conditions and soil types are shown to be in good agreement with Warrick’s semi-analytical solution and simulations using the HYDRUS-1D software package. Results obtained with the proposed method were far less dependent upon the grid spacing than the HYDRUS-1D finite element solutions.

Numerical solution of the Kirchhoff-transformed Richards equation for simulating variably saturated flow in heterogeneous layered porous media

A new numerical method was developed to accurately and efficiently compute a solution of the nonlinear Richards equation with a layered soil. In the proposed method, the Kirchhoff integral transformation was applied. However, in the Kirchhoff integral transformation approach, the transformed Kirchhoff head has dyadic characteristics at the material interface between different soil types. To avoid the dyadic characteristics at the material interface, a truncated Taylor series expansion was applied to the Kirchhoff head at the material interface and so the Kirchhoff head was replaced with a single pressure head value at the material interface. Accordingly, through the Taylor series expansion, a set of algebraic equations in the one-dimensional control volume finite difference discretized system formed a tridiagonal matrix system. Through a series of numerical experiments, the new method was compared to other numerical methods to determine its superiority. The results clearly demonstrated that the approach was not only more computationally efficient, but also more accurate and robust than other numerical methods. Computational performance was greatly enhanced with the proposed method, and which could be used to simulate complicated heterogeneous flow at a large-scale watershed or regional scale.

An improved iteration method for the numerical solution of groundwater flow in unsaturated soils

An improved Jacobi iterative method with two-side equilibration (TS-IJIM) is developed to model the process of unsaturated groundwater flow. The exponential model was adopted to linearize the nonlinear Richards equation (RE). Numerical discretization employing the finite difference method (FDM) produces a highly ill-conditioned system of linear equations. The TS-IJIM is proposed to reduce the matrix conditional number. Compared with two-side equilibration alone and the improved Jacobi iterative method alone, the proposed method effectively reduces the matrix conditional number, which is smallest when the resultant matrix has the same norm for each row. The TS-IJIM method is insensitive to the relaxation factor. Compared with the conjugate gradient method with a two-side equilibration algorithm, the TS-IJIM has better convergence in solving the highly ill-conditioned system. The conditional number of the ill-conditioned matrix can be reduced to 1.0 using the proposed TS-IJIM. Numerical results show that the proposed method can reduce the conditional number and improve convergence for solving unsaturated groundwater flow. Results also indicate that the proposed model is highly accurate and efficient in computing steady-state and transient 1D and 2D groundwater flow in unsaturated soils.

Reliability analysis of slope stability in earthen dams following rapid drawdown

This work was dedicated to assessing reliability of slope stability in earthen dams following rapid drawdown. Modeling of the problem physics was based on saturated-unsaturated seepage flow through the nonlinear Richards equation, under the assumption of two-dimensional approximation in steady state and transient state. The modeling of seepage taking place inside the earth dam infill materials included soil parameters, the negative pore water pressure-hydraulic conductivity relationship as well as the negative pore water pressure-volumetric water content. Slope stability was analyzed by using safety factor as evaluated by the Bishop's simplified method. Solution of the governing equations in terms of pore water pressure and fluid flow velocity was performed by using the modules SEEP/W and SLOPE/W of the Finite Element based software package GeoStudio. The drawdown phenomenon was then analyzed and its influence on dam stability assessed by using surface response method and Monte Carlo simulation. A design formula was proposed in order to bound the probability of failure associated to slope stability.

How does PTF Interpret Soil Heterogeneity? A Stochastic Approach Applied to a Case Study on Maize in Northern Italy

Soil water balance on a local scale is generally achieved by applying the classical nonlinear Richards equation that requires hydraulic properties, namely, water retention and hydraulic conductivity functions, to be known. Its application in agricultural systems on field or larger scales involves three major problems being solved, related to (i) the assessment of spatial variability of soil hydraulic properties, (ii) accounting for this spatial variability in modelling large-scale soil water flow, and (iii) measuring the effects of such variability on real field variables (e.g., soil water storage, biomass, etc.). To deal with the first issue, soil hydraulic characterization is frequently performed by using the so-called pedotransfer functions (PTFs), whose effectiveness in providing the actual information on spatial variability has been questioned. With regard to the second problem, the variability of hydraulic properties at the field scale has often been dealt with using a relatively simple approach of considering soils in the field as an ensemble of parallel and statistically independent tubes, assuming only vertical flow. This approach in dealing with spatial variability has been popular in the framework of a Monte Carlo technique. As for the last issue, remote sensing seems to be the only viable solution to verify the pattern of variability, going by several modelling outputs which have considered the soil spatial variability. Based on these premises, the goals of this work concerning the issues discussed above are the following: (1) analyzing the sensitivity of a Richards-based model to the measured variability of θ(h) and k(θ) parameters; (2) establishing the predictive capability of PTF in terms of a simple comparison with measured data; and (3) establishing the effectiveness of use of PTF by employing as data quality control an independent and spatially distributed estimation of the Above Ground Biomass (AGB). The study area of approximately 2000 hectares mainly devoted to maize forage cultivation is located in the Po plain (Lodi), in northern Italy. Sample sites throughout the study area were identified for hydropedological analysis (texture, bulk density, organic matter content, and other chemical properties on all the samples, and water retention curve and saturated hydraulic conductivity on a sub-set). Several pedotransfer functions were tested; the PTF‒Vereckeen proved to be the best one to derive hydraulic properties of the entire soil database. The Monte Carlo approach was used to analyze model sensitivity to two measured input parameters: the slope of water retention curve (n) and the saturated hydraulic conductivity (k0). The analysis showed sensitivity of the simulated process to the parameter n being significantly higher than to k0, although the former was much less variable. The PTFs showed a smoothing effect of the output variability, even though they were previously validated on a set of measured data. Interesting positive and significant correlations were found between the n parameter, from measured water retention curves, and the NDVI (Normalized Difference Vegetation Index), when using multi-temporal (2004–2018) high resolution remotely sensed data on maize cultivation. No correlation was detected when the n parameter derived from PTF was used. These results from our case study mainly suggest that: (i) despite the good performance of PTFs calculated via error indexes, their use in the simulation of hydrological processes should be carefully evaluated for real field-scale applications; and (ii) the NDVI index may be used successfully as a proxy to evaluate PTF reliability in the field.

Analytical solution of the one‐dimensional nonlinear Richards equation based on special hydraulic functions and the variational principle

SummaryA simple analytical solution of a one‐dimensional transient and nonlinear Richards equation is of great importance for estimating the hydraulic properties of soil and for precision irrigation. We developed simple approximate solutions in explicit and implicit forms for equations describing one‐dimensional, constant‐water‐head horizontal absorption and vertical infiltration for a special case of an exponential water‐retention curve and power‐law hydraulic conductivity, respectively, based on least‐action and variational principles. The profile of the soil water content (SWC) depended on the soil hydraulic parameter n associated with the pore‐size distribution, and the new solutions could be applied to most soils for 3.5 ≤ n ≤ 24.5. A comparison with numerical solutions calculated by an implicit‐difference scheme indicated that the approximate solution in explicit form accurately estimated the special hydraulic‐function parameters for calculating the SWC profile, cumulative infiltration, infiltration rate and wetting‐front distance in horizontal‐absorption experiments. An approximate solution in implicit form was also obtained with Nth order Taylor‐series expansions for an equation describing one‐dimensional, constant‐water‐head vertical infiltration. A large N was used to simulate SWC for soils with a large n, but this did not influence estimates of cumulative infiltration. The relation between infiltration rate and the inverse wetting‐front depth was not strictly linear for vertical infiltration, and the implicit method could be used to describe this nonlinear relation.Highlights A simple solution of the Richards equation is important to research soil water infiltration. Least‐action and variational principles were used to analyse soil water content distribution. Analytical solutions of a nonlinear Richards equation by special hydraulic functions were proposed. A method was proposed to predict the infiltration and estimate soil hydraulic parameters.

Infiltration Model Based on Traveling Characteristics of Wetting Front

Core Ideas A new infiltration model based on characteristics of wetting front was developed. The model describes evolution of water content by ordinary differential equation. Infiltration rate is used as input, and wetting front profile and time are back‐solved. The model captures the evolution of wetting front and infiltration amount. This study presents a one‐dimensional infiltration model based on the traveling characteristics of wetting front (TCWF) model. The TCWF model predicts the infiltration process by solving an ordinary differential equation linking infiltration rate and wetting front profile, which avoids numerically solving the highly nonlinear Richards equation and eliminates associated numerical difficulties. This is achieved by an assumption of normalized equality between flux and soil moisture within wetting front before the time when infiltration rate reaches steady value and traveling wave property of wetting front in the equilibrium state. With the two traveling characteristics, the TCWF model uses infiltration rate as input and back solves wetting front profile, time, and cumulative infiltration. We examine the accuracy and capability of TCWF model by comparing 48 infiltration cases with different setups of soil profile and infiltration conditions to the numerical solutions of the Richards equation for vertical infiltration and exact solutions for horizontal infiltration. Results demonstrate that the TCWF model can effectively capture the evolution of wetting front and predict infiltration. The TCWF model can be applied to any soil type, soil hydraulic model, and infiltration on any slope angle. Good agreements with maximum percent bias for cumulative infiltration of 3.525% are found in the vertical infiltration cases with initially dry and initially wet conditions and in the horizontal case with initially dry conditions. In the absence of gravity, the TCWF model slightly overestimates horizontal infiltration, especially under initially wet conditions. The developed TCWF model provides a simple and accurate approach to simulate infiltration.

NUMERICAL MODELING OF UNSATURATED LAYERED SOIL FOR RAINFALL-INDUCED SHALLOW LANDSLIDES

In this paper, a pioneer study on numerical modeling of rainfall-induced shallow landslides in unsaturated layered soil using the variably saturated flow equation is presented. To model the shallow landslides, the infinite slope stability analysis coupled with the hydrological model with the consideration of the fluctuation of time-dependent pore water pressure and Gardner equation for soil water characteristic curve was developed. A linearization process for the nonlinear Richards equation to deal with groundwater flow in unsaturated layered soil is derived using the Gardner model. To solve one-dimensional flow in the unsaturated zone of layered soil profiles, flux conservation and the continuity of pressure potential at the interface between two consecutive layers are considered in the numerical discretization of the finite difference method. The validity of the proposed model is established in three numerical problems by comparing the results with the analytical and other numerical solutions. Application examples have also been conducted. Obtained results demonstrate that the fluctuation of pore water pressure in unsaturated layered soil dominates slope stability of landslides and the lowest factor of safety may occur at the interface between two consecutive layers. The findings observed in this study are a fundamental contribution to environmental protection engineering for landslides in areas with higher occurrence and vulnerability to extreme precipitation.

Transient Modeling of Flow in Unsaturated Soils Using a Novel Collocation Meshless Method

In this paper, a novel meshless method for the transient modeling of subsurface flow in unsaturated soils was developed. A linearization process for the nonlinear Richards equation using the Gardner exponential model to analyze the transient flow in the unsaturated zone was adopted. For the transient modeling, we proposed a pioneering work using the collocation Trefftz method and utilized the coordinate system in Minkowski spacetime instead of that in the original Euclidean space. The initial value problem for transient modeling of subsurface flow in unsaturated soils can then be transformed into the inverse boundary value problem. A numerical solution obtained in the spacetime coordinate system was approximated by superpositioning Trefftz basis functions satisfying the governing equation for boundary collocation points on partial problem domain boundary in the spacetime coordinate system. As a result, the transient problems can be solved without using the traditional time-marching scheme. The validity of the proposed method is established for several test problems. Numerical results demonstrate that the proposed method is highly accurate and computationally efficient. The results also reveal that it has great numerical stability for the transient modeling of subsurface flow in unsaturated soils.

Analytical models of infiltration and redistribution for unsaturated flow in soils with vertically non-uniform saturated hydraulic conductivity

The Richards equation describes unsaturated water movement in soils. However, analytical solutions to the nonlinear Richards equation are possible only for simple and idealized conditions. In this paper, sharp front (SF) models are developed for vertically non-uniform soils, with the saturated hydraulic conductivity continuously decreasing with depth according to a power-law function. Analytical expressions for evolutions of water content profiles and fluxes past a control plane are developed for the three surface boundary conditions of constant saturation, applied surface flux, and redistribution that are applicable during and after rainfall events. Comparisons of model results with infiltrometer experimental data from Purdue University are provided for the special case of falling head conditions at the soil surface. Subsequently, model results are compared with numerical solutions of the Richards equation for different levels of vertical heterogeneity for rainfall events of finite duration so that all the three boundary conditions are employed. Results show that soil moisture fronts are steeper with increasing non-uniformity, making sharp front approximations more accurate for pre- and post-ponding conditions, but less so for redistribution. The developed model results are found to provide reasonable approximations to numerical solutions of the Richards equation for vertically non-uniform soils.

On the construction of a functional solution method for the infiltration in porous media problems

In this paper, we consider a transient vertical one-dimensional flow problem of water in unsaturated porous media, modelled by the non-linear Richards equation. Constitutive relations of Van Genuchten will be employed and Padé approximants are used to represent the hydraulic capacity and conduction in a simplified fashion. From the proposed methodology a construction of a functional solution is presented with the objective to define an initialisation of a recursive scheme in the spirit of Adomian decomposition. This solution is optimised and evaluated using the governing equation for a self-consistency test. The results are presented for some soil types and its related soil parameters, that are reported in the literature.

Numerical Solutions for Groundwater Flow in Unsaturated Layered Soil with Extreme Physical Property Contrasts

Abstract In this paper, the numerical solutions for groundwater flow in unsaturated layered soil using the Richards equation are presented. A linearisation process for the nonlinear Richards equation to deal with groundwater flow in unsaturated layered soil is derived. To solve one-dimensional flow in the unsaturated zone of layered soil profiles, flux conservation and the continuity of pressure potential at the interface between two consecutive layers are considered in the numerical model. In addition, a novel method, named the dynamical Jacobian-inverse free method, incorporated with a two-side equilibration algorithm for solving ill-conditioned systems with extreme contrasts in hydraulic conductivity is proposed. The validity of the model is established in numerous test problems by comparing the numerical results with the analytical solutions. The results show that the proposed method can improve convergence and numerical stability for solving groundwater flow in unsaturated layered soil with extreme contrasts in hydraulic conductivity.

Solving the nonlinear Richards equation model with adaptive domain decomposition

Modeling the transport processes in a vadose zone plays an important role in predicting the reactions of soil biotopes to anthropogenic activity, e.g. modeling contaminant transport, the effect of the soil water regime on changes in soil structure and composition, etc. Water flow is governed by the Richards equation, while the constitutive laws are typically supplied by the van Genuchten model, which can be understood as a pore size distribution function. Certain materials with dominantly uniform pore sizes (e.g. coarse-grained materials) can exhibit ranges of constitutive function values within several orders of magnitude, possibly beyond the length of real numbers that computers can handle. Thus a numerical approximation of the Richards equation often requires the solution of systems of equations that cannot be solved on computer arithmetics. An appropriate domain decomposition into subdomains that cover only a limited range of constitutive function values, and that will change adaptively, reflecting the time progress of the model, will enable an effective and reliable solution of this problem. Parts of this problem have already been treated in our work Kuraz et al. (2014). This paper focuses on improving the performance of a nonlinear solver by locating the areas with abrupt changes in the solution.