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

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 lumping, time weighting, conductivity and moisture content characteristics, boundary conditions, seepage faces, and the extent of fully saturated zones in the soil. Previous strategies for enhancing the performance of the Picard 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 Jacobian matrix, and it is suggested that both methods can be effectively implemented and used in numerical models of Richards' equation.