A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems

Abstract

Picard iteration is a widely used procedure for solving the nonlinear equation governing flow in variably saturated porous media. The method is simple to code and computationally cheap, but has been known to fail or converge slowly. The Newton method is more complex and expensive (on a per-iteration basis) than Picard, and as such has not received very much attention. Its robustness and higher rate of convergence, however, make it an attractive alternative to the Picard method, particularly for strongly nonlinear problems. In this paper the Picard and Newton schemes are implemented and compared in one-, two-, and three-dimensional finite element simulations involving both steady state and transient flow. The eight test cases presented highlight different aspects of the performance of the two iterative methods and the different factors that can affect their convergence and efficiency, including problem size, spatial and temporal discretization, initial solution estimates, convergence error norm, mass lureping, time weighting, conductivity and moisture content characteristics, boundary conditions, seepage faces, and the extent of fully satarated zones in the soft. Previous strategies for enhancing the performance of the heard and Newton schemes are revisited, and new ones are suggested. The strategies include chord slope approximations for the derivatives of the characteristic equations, relaxing convergence requirements along seepage faces, dynamic time step control, nonlinear relaxation, and a mixed Picard-Newton approach. The tests show that the Picard or relaxed Picard schemes are often adequate for solving Richards' equation, but that in cases where these fail to converge or converge slowly, the Newton method should be used. The mixed Picard-Newton approach can effectively overcome the Newton scheme's sensitivity to initial solution estimates, while comparatively poor performance is reported for the various chord slope approximations. Finally, given the reliability and efficiency of current conjugate gradient-like methods for solving linear nonsymmetric systems, the only real drawback of using Newton rather than Picard iteration is the algebraic complexity and computational cost of assembling the derivative terms of the JacobJan matrix, and it is suggested that both methods can be effectively implemented and used in numerical models of Richards' equation. great popularity because it is the most intuitive linearization of Richards' equation, is computationally inexpensive on a per-iteration basis, and preserves symmetry of the discrete system of equations. However, the method may diverge under certain conditions, as has been observed empirically [e.g., Huyakorn et al., 1984; Celia et al., 1990] and verified theoretically [Aldama and Paniconi, 1992]. The Newton scheme, also known as Newton-Raphson iteration, yields nonsymmetric system matrices and is more complex and expensive than Picard linearization, though it achieves a higher rate of convergence and can be more robust than Picard for certain types of proble...

Incipient subsurface heterogeneity and its effect on overland flow generation – insight from a modeling study of the first experiment at the Biosphere 2 Landscape Evolution Observatory

“…The subsurface flow module in CATHY solves the 3-D Richards equation describing flow in variably saturated porous media (Paniconi and Putti, 1994), while the surface flow module solves the diffusion wave equation describing surface flow propagation over hillslopes and in stream channels identified using terrain topography and the hydraulic geometry concept (Orlandini and Rosso, 1998). Surface-subsurface coupling is based on a boundary condition switching procedure that automatically partitions potential fluxes (rainfall and evapotranspiration) into actual fluxes across the land surface and calculates changes in surface storage.…”

Assessment of the Impact of Subsurface Agricultural Drainage on Soil Water Storage and Flows of a Small Watershed

“…In their sensitivity analysis study of CATHY to the soil hydrodynamic properties, Muma et al [53] noticed that the saturated hydraulic conductivity of the deeper layers (fifth group of layers) had a significant impact on to drain discharge and outlet of the micro-watershed flow. Furthermore, they revealed that the vertical saturated hydraulic conductivity in the two surface layers (first group of layers) as well as the vertical and lateral saturated hydraulic conductivity in the layers where the subsurface drains are located deserved special attention due to their strong interaction with other parameters with regards to drain discharge.…”

“…Isotropic and anisotropic properties of soil saturated hydraulic conductivity (*) analysed. In their sensitivity analysis study of CATHY to the soil hydrodynamic properties, Muma et al [53] noticed that the saturated hydraulic conductivity of the deeper layers (fifth group of layers) had Notes: * All hydraulic conductivity values are expressed in m/s. In each column, the value on the left of the slash is the horizontal hydraulic conductivity (X and Y), while the value on the right is the vertical hydraulic conductivity (Z).…”

“…Since the storage and conductivity terms strongly depend on pressure, Equation (2) is mostly nonlinear and, as a result, is linearized by Picard's or Newton's iterative methods [53]. Thus, the partial differential equation, which mathematically describes the flow in porous media, is:…”

Efficient method for simulating gravity‐dominated water flow in unsaturated soils

“…Modeling of water flow in unsaturated porous media is usu-, ally based on Richards' equation [Richards, 1931], and with a few exceptions, existing numerical models use Picard or Newton iteration [Paniconi and Putti, 1994]. Iteration-to-iteration oscillation often arises and causes slow convergence or failure to converge, especially when simulating large and rapidly changing infiltration flows .…”