A Functional Iteration Technique for Solving the Richards Equation Applied to Two‐Dimensional Infiltration Problems

Abstract

In the solution of nonlinear parabolic partial differentiM equations, such as the Richards equation, classical implicit schemes often oscillate and fail to converge. A fully implicit scheme has been developed along with a functional iteration method for solving the system of nonlinear difference equations. Newton's iteration technique is mathematically the most preferable of all functional iteration methods because of its quadratic convergence. The Richards equation, Newton-linearized with respect to relative permeability and saturation as functions of capillary pressure, is particularly aided by this new approach for problems in which saturations vary rapidly with time (infiltration fronts, cone of depression near a well bore, and so forth). Although the computing time is almost twice as long for a time step with Newton's iteration scheme, the smaller time truncation than that of classical implicit schemes and the stability in cases in which classical schemes are unstable permit the use of much larger time steps. To demonstrate the method, heterogeneous (layered) soil systems are used to simulate sharp infiltration fronts caused by ponding at the soil surface. When water filters into or through a soil, it replaces or is replaced by air. It may often bepossible to ignore the flow of the air phase because of the relatively small resistance to the flow of air. This assumption is made in the present paper so that the analysis of fluid flow in a soil saturated partially with water and partially with air reduces to one-phase flow, i.e., flow of water only. In such cases it is convenient to consider the pressure of the air to be constant everywhere and usually atmospheric [Corey, 1969]. Under these conditions the Richards equation [Richards, 1931], a combination of continuity and Darcy's law for the water phase only, applies. The present paper proposes a numerical scheme for a two-dimensional x, z or r, z form of the Richards equation (equation 1) suitable for infiltration problems in heterogeneous (layered) and anisotropic soils. A comparison of numerical and experimental results will be made in a subsequent paper along with a demonstration of solutions to practical problems. To date one-dimensional numerical results of the Richards equation are numerous and have been reviewed by Freeze [1969]. Two-dimensionM applications are less numerous. Most of them relate to 'free surface' flow problems integrating the saturated and unsaturated zones into a hydrodynamic unity [Rubin, 1968; Hornberger et al., 1969; Taylor and Luthin, 1969; Cooley, 1970; Verma and Brutsaert, 1970]. A unique threedimensional model was constructed (R. A. Freeze, unpublished report, 1970) for the treatment of saturated-unsaturated transient flow in small nonhomogeneous anisotropic geologic basins. In all of the preceding studies except that of Freeze (unpublished report, 1970) homogeneous soils are used in the analysis. All two-dimensional studies are restricted to cases in which saturations vary rather smoothly in space and time; i.e., sharp infiltra...

Algorithms for solving Richards' equation for variably saturated soils

“…The iterative alternating direction implicit, the line successive over relaxation, and the strongly implicit method algorithms have been used with varying degrees of success [Freeze, 1971a]. Brutsaert [1971] was one of the first authors to present a finite difference algorithm which combined the mixed form of Richards' equation with Newton iteration to effectively deal with steep wetting fronts. Recently, a mixed form of Richards' equation using modified Picard iteration and a preconditioned conjugate gradient solver has shown much promise in modeling unsaturated problems that contain steep wetting fronts [Bouloutas, 1989;Celia et al, 1990].…”

Physically based modeling in catchment hydrology at 50: Survey and outlook

“…Accuracy and mass conservation can be improved by using alternative formulations of RE or by applying appropriate primary variable transformations [ Celia et al ., ; Kirkland et al ., ; Rathfelder and Abriola , ; Pan and Wierenga , ; Diersch and Perrochet , ; Williams et al ., ]. Solving RE requires that the equation be linearized, and Newton‐based iterative schemes, including less accurate approximations such as the Picard method, are commonly used for this purpose [ Brutsaert , ; Huyakorn et al ., ; Paniconi and Putti , ; Gustafsson and Söderlind , ; Casulli and Zanolli , ; Lott et al ., ], although noniterative approaches have also been proposed [ Paniconi et al ., ; Kavetski et al ., ; Ross , ; Crevoisier et al ., ]. The convergence behavior of an iterative scheme at each time step of the RE solution can be advantageously used to adapt the step size during the simulation, while noniterative schemes can rely on local truncation error estimates to dynamically control the step size [ D'Haese et al ., ].…”

Self‐adaptive hierarchic finite element solution of the one‐dimensional unsaturated flow equation

“…If the nodeless variables are ordered so that they follow all the nodal values of h, the matrix equation (8) can be rewritten after enrichment with quadratic basis functions as In equation(l2) the matrix M I , is a tridiagonal matrix composed of the same terms as the equivalent Galerkin formulation using only linear basis functions. To formulate M , , and M, the corresponding matrices A and B need to be examined.…”