Nonlinear Effects on the Convergence of Picard and Newton Iteration Methods in the Numerical Solution of One-Dimensional Variably Saturated – Unsaturated Flow Problems

Finite element discretization of the pressure head form of the Richards equation leads to a nonlinear model, which yields numerical convergence difficulties. When the numerical solution to this problem has either an extremely sharp moving front, infiltration into dry soil, flow domains containing materials with spatially varying properties, or involves time-dependent boundary conditions, the corrector iteration used in many time integrators can terminate prematurely, which leads to incorrect results. While the Picard and Newton iteration methods can solve this problem through tightening the tolerances provided to the solvers, there is a more efficient approach to overcome the convergence difficulties. Four tests examples are examined, and each test case is solved with five sufficiently small tolerances to demonstrate the effectiveness of convergence. The numerical results illustrate that the methods greatly improve the convergence and stability. Test experiments show that the Newton method is more complex and expensive on a per iteration basis than the Picard method for simulating variably saturated–unsaturated flow in one spatial dimension. Consequently, it is suggested that the resulting local and global mass balance is exact within the minimum specified accuracy.


Introduction
Fluid flow through variably saturated-unsaturated porous media is usually described by the classical Richards equation [1].It is defined by coupling a statement of mass conservation with Darcy's equation.It is a highly nonlinear partial differential equation that can be represented in three standard forms, depending on whether pressure, moisture, or both are used as dependent variables.For one-dimensional vertical flow, the pressure head-based Richards equation is written as where, Ψ is the pressure head (L), t is time (T), z denotes the vertical distance from reference elevation, which is assumed positive upward (L), K(Ψ) is the hydraulic conductivity (LT −1 ), C(Ψ) = dθ dΨ is the specific fluid capacity (L −1 ), and θ is the moisture content.
Numerical schemes for solving Equation (1) on the basis of grid resolution techniques, and/or time discretization procedures including other dynamic approaches for different forms of Richards equation are available.The choices of these techniques greatly influence the computational efficiency, numerical stability, and accuracy of water flow problems.The moisture content of the Richards equation is limited to unsaturated flow problems due to θ not being a smooth function for layered soils, and discontinuous near saturation.As a result, severe numerical convergence difficulties, numerical oscillations, and stability problems result in the failure of subsurface solvers.Thus, this form may be useful only for homogeneous flow conditions.On the other hand, a pressure head formulation of Richards equation is continuous in both saturated and unsaturated zones, and can be used for homogeneous and non-homogeneous soils.However, numerical approaches based on this form can suffer large mass balance errors [2][3][4].Small step size coupled with mass lumping effectively ensures the improvement of the mass balance of the pressure head-based Richards equation [2].
Due to the shape of the soil water retention and hydraulic conductivity curves, the Richards equation is highly nonlinear.For the less sharp region of the moisture content or moisture capacity, many iterations are required to achieve a robust and accurate solution of the Richards equation with relatively large tolerance.For the portion of the soil water curve where the profile is very sharp, a much smaller specified accuracy should be used to avoid the numerical errors as well as computational challenge [5].Therefore, the numerical solution of the Richards equation is very sensitive to the value of the convergence tolerance.An adaptive and heuristics time-stepping scheme [6], a mass-conservative algorithm [2], and a nested-Newton algorithm [7] with small nonlinear tolerance have been successfully applied for different soil properties with different initial and boundary conditions depending on the formulation of the Richards equation.Besides a mining geostatical approach [8] for disordered soil's structure, a fully analytical model in case of nonstationary of unsaturated flow variables has been introduced [9].Therefore, the flow problems that involve extremely sharp front in space-and/or time-varying boundary conditions flow into layered soil, making them difficult to solve accurately unless tight tolerance is used.
The highly nonlinear model ( 1) is solved using common constitutive relations and standard approximation methods such as finite differences [2,3,10], finite elements [2,[11][12][13], and finite volumes [12,14,15].Standard algorithms that employ the numerical solution of sharp front problems with poor accuracy are computationally expensive and lead to less accurate results.The solution of the resulting large nonlinear algebraic systems that arise in implicit numerical discretizations of the Richards equation has been the subject of significant research, and hence, one needs efficient and robust linearization techniques that maintain not only the accuracy of the solution, but also its mass conservation property.Typical linearization methods, such as the Picard, Newton, fast-secant and relaxation methods, as well as non-iterative methods, such as the implicit factored schemes and three-level Lees schemes, have been proposed [4,11,[16][17][18].
In the process of linearization, two iterative schemes, Picard and Newton, are characterized by different convergence behavior.The Picard iterative scheme is more popular than the Newton due to its simplicity and satisfactory performance [2].Theoretically, the Newton method converges one order faster than the Picard method, but under certain flow conditions, the Picard method is more efficient than the Newton method [16].Generally, the Picard method is less expensive per iteration basis, and preserves the symmetry of the discrete system of equations.However, it suffers from slow convergence, or may diverge for cases of gravity drainage, complex time-varying boundary conditions, strongly nonlinear characteristic equations, and saturated-unsaturated interfaces [16].The Newton method has better convergence properties, but creates nonsymmetric system matrices and involves the computation of derivatives, which leads to a complex and expensive linearization technique.This study concerns the comprehensive behavior of two linearization methods, Picard and Newton, with tighter accuracy for efficiently solving the Richards equation in saturated-unsaturated porous media based on backward Euler finite difference in time and finite elements in space with mass lumping.
The objectives of this work are: to assess the lack of robustness of solution of the Richards equation when using common solution approaches with very small tolerances, to investigate the effect of different tolerances, and to compare a set of alternative tolerances for a wide range of porous media conditions to test robustness and efficiency with minimum accuracy limit.

Finite Element Model
To develop a finite element approximation of (1), a Galerkin finite element method is used to discretize the spatial domain with linear basis function.Solutions are obtained by solving the discretized Richards equation using a mass-lumping linear Galerkin finite element method, which has been proved to be a mass conservation formulation [2].The time discretization is carried out with implicit backward Euler finite difference method, and mass lumping is used for the mass accumulation term.Discretization yields a system of nonlinear equations: where Ψ is the vector of the pressure head at each node, k denotes the time step, and A, F, b and q are the stiffness matrix, mass matrix, gravitational gradient component, and Darcy's flux boundary conditions, respectively.

Linearization Techniques
To obtain a cost-effective solution of the system of nonlinear Equation (2), Picard iterative and a more strongly convergent Newton iterative schemes are employed.The later one is designed especially to avoid the convergence difficulties usually encountered when certain nodes are being desaturated.However, this technique is computationally more expensive and algebraically complex than the Picard scheme.

Newton and Picard Schemes
Letting superscript (m) be an iteration level, using (2), the Newton scheme [16] can be written as where h = Ψ k+1, (m+1) − Ψ k+1, (m) , and the Jacobian for the system is The simplest formulation of the Picard scheme [16] of the Richards equation is

Characteristic Equations
Solving Richards equation requires constitutive relationships relating the dependent variable pressure head and the nonlinear terms, such as moisture content, moisture capacity, and conductivity.The most commonly used soil characteristic models, van Genuchten [19] and Brooks-Corey [20], are employed in this work, and these models are illustrated below.

Van Genuchten Model
The most frequently used constitutive relationships were derived by van Genuchten [19] and are given by: where, m = 1 − 1 n is a pore-size distribution index, θ s is the saturated moisture content, and θ r is the residual moisture content.

The Brooks-Corey Model
Another common form of the constitutive relationships proposed by Brooks and Corey [20] are given by: where ψ d = − 1 α is the bubbling or air entry pressure head (L).

Implementation
Numerical results were obtained by implementing the Picard and Newton iteration techniques into the computer code CATHY (CATchment Hydrology) [21].It is a coupled physically-based spatially-distributed model for surface-subsurface simulations.The model simulates the three-dimensional equation for subsurface flow in variably saturated porous media (i.e., Richards equation) using Galerkin finite elements in space, a weighted finite difference scheme in time, and linearization via Newton or Picard iteration.In the Picard and Newton iterative schemes, convergence is achieved when the termination criterion Ψ k+1, (m+1) − Ψ k+1, (m) ≤ Tol is satisfied, where Tol is the specified tolerance/accuracy, and in this work, five sufficiently small enough values (10 −2 , 10 −3 , 10 −4 , 10 −5 and 10 −6 ) are used.In all of the runs, the time-step size is adaptively adjusted according to the convergence behavior of the previous nonlinear iteration [16] in order to optimize convergence and CPU efficiency.The simulation begins with a time step size of ∆t 0 , and continues until reach the ending time T max .At each time step of the simulation, the current time-step size is increased by a factor of ∆t mag (=1.20) to the specified maximum size of ∆t max if convergence at the previous iteration is achieved in fewer iteration than maxit 1 (=5), but it is left unchanged if convergence requires between maxit 1 and maxit 2 (=8).When the scheme is not convergent within the iteration maxit 2 , the time step is decreased by a reduction factor of ∆t red (=0.5) to a minimum of ∆t min .If convergence is not achieved (maxit exceded), the time-step size is repeated ('back-stepping') using a reduced time-step size (factor ∆t red , to a minimum of ∆t min ).The values of the ∆ts and maxits are chosen.The main advantage of the dynamic time marching is that it avoids undesirable abrupt changes in time step.Linear solvers for large sparse nonsymmetric systems arising from the Newton scheme are generally inefficient, but presently consistent, and efficient conjugate gradient-type algorithms are available for solving such systems.Biconjugate gradient stabilized algorithm (BICGSTAB), minimum residual algorithm (GRAMRB), generalized conjugate residual method (GCRK), and transpose-free quasi-minimal residual algorithm (TFQMR) are reliable, large, sparse nonsysmmetric solvers.Detailed descriptions of the various algorithms can be found in the works Axelsson, Pini et al., van der Vorst, and Freund [22][23][24][25].The solutions of the linear systems obtained by using BICGSTAB algorithm in modeling flow problems have shown best performance [17].Another important feature of this solution's employed algorithms is the use of an analytical differentiation for the evaluation of the gradient of the soil moisture characteristic curves.To measure the accuracy of the solution of the subsurface solver, mass balance error (MBE) is calculated.According to the mass continuity equation, the MBE is defined as the difference between the net amount of water added to the system, and the change in the amount of water stored in the system after a given elapsed time, and it is evaluated by the following formula [2].
with N = E + 1 nodes {z 0 , z 1 , z 2 , . . . ,z E }, constant nodal spacing ∆z is considered, and q o and q N are boundary fluxes calculated from the finite element equations associated with the boundary nodes z 0 and z N .For this study, the simple procedure given above has been proved to be effective and sufficient.All of the numerical simulations were executed on personal machine, Dell INSPIRON 2.5 GHz with 4 GB RAM of intel Core 2 processor.

Numerical Results
Four one-dimensional test simulations with different soil characteristics were performed for evaluating the performance of two iterative schemes for five sufficiently small accuracies.The numerical results needed to meet the objective of this work; we focus our attention on the analysis of how much the small tolerance will affect the efficiency, accuracy, and convergence properties of the numerical solution to the Richards equation.However, the performance of the simulations is more difficult to judge, as there are a number of factors that control the efficiency and robustness of the model.Such factors include total central processing unit (CPU), total number of time steps, average nonlinear (NL) iterations (Iter) per time step, cumulative mass balance error (Cum.MBE), the number of backstepping, etc.The most significant parameters, depending on the physical and numerical characteristics of the specific simulations, will be selected to better convey the performance of the Picard and Newton iterative schemes.

Problem 1
This one-dimensional test problem consists a soil column of 0.3 m length discretized with a fine vertical resolution ∆z = 0.0024 m.This test case is a standard test problem [26], and substantially easier than remaining test problems because the domain is much shorter, media is nonuniform, and saturated conditions are not developed.However, it is a common test case and very helpful in clearly illustrating the analysis of nonlinear performance of the proposed methodology.The initial pressure head is −10 m, and constant pressure head boundary conditions are applied at the bottom and top of the soil column; these are −10 m and −0.75 m, respectively.The soil hydraulic properties are described by the van Genuchten model with θ s = 0.368, θ r = 0.102, α = 3.35/m, n = 2.0, and K s = 7.970 m/day.As a result, simulation conditions yield a difficult sharp-front problem.
Computed pressure head profiles for different tolerances of two iterative schemes are presented in Figure 1.Both the constant pressure head and initial condition of this test problem are easily noted.Figure 1 shows that it is very difficult to distinguish the solutions obtained with various tolerances for two iterative schemes, and all of the profiles are very similar to the published results [26][27][28][29].However, slight numerical oscillations are observed in the bottom of the soil column for the Picard scheme with error tolerance 10 −2 m.The computational performance of this test problem of two linearization techniques can be compared with the results of different tolerances, as shown in Table 1.It is clear that the cumulative mass balance errors for all of the runs are very close and remarkably small throughout the entire simulation, confirming the consistency of the schemes and reliability of the accuracies.The measures of other computational efforts used in this study are the total number of time steps and average nonlinear iterations per time step.Usually, these results vary with the size of time steps, because the step sizes are strongly influenced by the convergence of the nonlinear solver.When we reduce the tolerance by a factor of 10 −1 m, the number of time steps increases approximately 50-80%.Central processing unit (CPU) time is a common measure, but it is dependent on machine, compiler, program, system load, and algorithm.It is the most meaningful measure for evaluating the relative efficiency for solving the Richards equation, and the results in Table 1 indicate that CPU time is greatly influenced by the reduction of tolerance.For this simple test, the Picard method is characterized by a lower (approximately one-fifth for the same tolerance) computational cost with respect to the Newton technique because of the differences in the per iteration costs of the two approaches, as can be seen in the last row of Table 1.Almost the same number of average nonlinear iterations is required for all of the cases.We note that both the schemes resulted in significantly fewer backstepping occurrences.Convergence behavior and MBE plots for the Picard and Newton schemes are presented in Figures 2 and 3, respectively.The results of the two iterative schemes with tolerance 10 −3 m are indistinguishable from all of the solutions within the simulation and hence, we can say that it has sufficient error tolerance to achieve an accurate, efficient, and robust solution under such flow conditions.

Problem 2
This one-dimensional test case concerns a difficult vertical infiltration problem of a 10 m-high soil column.It will be solved, and has already been analyzed in detail by many researchers [26,[28][29][30].This soil column is parameterized by the van Genuchten model with the soil parameters θ s = 0.301, θ r = 0.093, α = 5.47/m, n = 4.264, and K s = 5.040 m/days.
Constant head boundary conditions at both the top (ψ(10, t) = 0.1) and bottom (ψ(0, t) = 0.0) boundaries, and a hydrostatic equilibrium initial condition (ψ(z, 0) = −z), are applied.These forcing conditions, along with the constitutive relationships, lead to the development of a sharp infiltration front and induce large gradients in the solution.This type of problem provides a rigorous test case for time integrators, and is well suited for the analysis of numerical convergence and efficiency with accuracy.As a result, it is a meaningful test problem to measure the limit of accuracy.
Due to the expense involved in accurate computing, a sufficient number of nodes in a grid is important to approximate the Richards equation.Since our objective is to achieve accurate results with the highest possible tolerance, the dense grid solutions for this test problem were made on a uniform grid of 401 nodes.
The solutions from these dense grid simulations at time T = 17,280 s for the Picard and Newton methods with different tolerances are shown in Figure 4.A sharp infiltration front between the saturated and unsaturated zones and the drained-to-equilibrium initial condition are noted, which are the trademark of this test problem.It is evident from these figures that the same accuracy is obtained for each of the five tolerances, and the results show good agreement with that of the published reports [26,[28][29][30].The results obtained by two iterative techniques with different tolerances are compared in Table 2 and Figures 5 and 6.A low mass balance error is necessary, but is not a sufficient condition to ensure the accuracy of the solution.In other words, accurate solutions ensure a small mass balance error.The MBE results for the two iterative schemes with five tolerances confirm the accurate numerical solutions.Errors are not significantly decreased for all of the runs if the tolerances are decreased, which are almost same for all of the runs performed in this work.For the Picard and Newton runs, the reduction in the magnitude of the tolerance in the nonlinear solver does not significantly affect the overall accuracy of the solution.There are very little differences for the nonlinear iterations per time step.For both the Picard and Newton schemes, approximately 1.5 times and three times more steps are required to converge to the final solution compared with each of the previous head tolerances, respectively.The same inclinations are observed in CPU time for all of the given levels of iteration tolerances, which drastically increases for the convergence of the Newton scheme.The oscillations of nonlinear convergence (Figures 5 and 6) and the number of backsteppings are comparatively very high for the case of 10 −2 m, which can be attributed to a sharp infiltration front produced at top of the soil column.Therefore, as far as the accuracy of the solution of the Richards equation is concerned, 10 −3 m is enough as a level of tolerance to achieve a reliable and robust solution.

Problem 3
In order to attain our objective, our present test study consists of one-dimensional vertical infiltration with redistribution in a 5 m soil column [29][30][31] Comparison of pressure head profiles in the Picard and Newton methods by using all values of the head tolerances are shown in Figure 7. Due to the dynamic boundary conditions at the top of the domain, there is rapid infiltration of water from the surface, followed by a period of redistribution of the water.The simulated results in each of the two schemes with five tolerances closely agree with the previous studies [29][30][31].A major problem in solving the Richards equation is to maintain good MBE concerning its nonlinear nature, especially for flow involving extremely sharp two-front and rapid infiltration of water from the surface with dynamic boundary conditions.A small mass balance does not guarantee an appropriate numerical solution; for an accurate solution, it is usually required to satisfy the mass conservation property, and the evaluated MBE results presented in Table 3 attest to this proposition.It can be seen from the computational results for all of the runs that the two iteration methods precisely conserved the mass.Very closed MBE results are obtained for all cases; however, some differences are observed in terms of average nonlinear iterations per time steps, except the case of 10 −2 m, it should be noted from Table 3 that both numerical techniques have taken the largest average time step size during the simulation; as a result, the number of iterations becomes smallest for the head tolerance 10 −3 m.For this test problem, the Picard and Newton schemes need on average four to five iterations per time step throughout the simulation period for all of the head tolerances.This explains why remarkable differences are not shown by the Picard and Newton schemes in terms of convergence behavior, as we observed on the plots of Figures 8 and 9 (including MBE behavior).The number of backstepping increased dramatically compared with the previous examples due to the strong nonlinear nature of the soil hydraulic properties with time-varying surface boundary conditions.Also, both the Picard and Newton iterative methods for the case of 10 −2 m failed to converge several times (56 for Picard and 48 for Newton) during the course of simulation.It can be inferred from these discussions that the head tolerance of 10 −3 m bears consistent head distribution profiles, as well as perfect mass conservation nature.

Problem 4
This case involves vertical drainage through a layered soil from initially saturated conditions.At time t = 0, the pressure head at the base of the column is reduced from 2 m to 0 m.During the subsequent drainage, a no-flow boundary condition is applied to the top of the column.This problem is considered to be a challenging test for numerical methods, because a sharp discontinuity in the moisture content occurs at the interface between two material layers [7,15,32].
During downward draining, the middle coarse soil tends to restrict drainage from the upper fine soil, and high saturation levels are maintained in the upper fine soil for a considerable period of time.The Brooks-Corey model is used to prescribe the pressure-moisture relationship.The hydraulic properties of the soils are given in Table 4.The soil profile is Soil 1 for 0 < z < 60 cm and 120 cm < z < 200 cm, and Soil 2 for 60 cm < z < 120 cm.A Dirichlet boundary condition is imposed at the base of the bottom boundary.Simulations are performed by the Picard and Newton methods on a fine mesh of 150 elements with five very small head tolerances.The evolution of saturation distribution in the soil column while drainage occurs is illustrated in Figure 10 for all of the cases of head tolerances at an elapsed time of 1,050,000 s (approximately 12 days).All of the solutions are in excellent agreement with the published results [7,15,32], except the Newton runs for error tolerance 10 −2 m.Desaturation takes place at early times in both the middle medium sand and upper fine sand layers, as clearly shown in this plot.A physical implication of the Richards equation, as this figure demonstrates, is that air is freely available to all parts of the water flow system, irrespective of the saturation distribution.At the intermediate times, desaturation in the upper fine sand layer decelerates significantly, as the middle medium sand keeps desaturating.It is clear from the figure that the base of the middle medium sand is mostly desaturated, and the drying front has begun to move downward into the lower fine sand layer.At later times, desaturation occurs very slowly in the upper fine sand layer.Finally, a new hydrostatic saturation distribution will be formed in the column.A compelling feature of this simulation is that the middle medium sand layer tends to limit drainage from the overlying fine sand.High saturations are presented in the upper fine sand layer for significant periods of time.This test example demonstrates the efficacy of coarser-grained layers as capillary barriers in unsaturated flow.
Numerical simulations for layered soil under vertical drainage are generally much more complicated than those under homogeneous porous media, and require more CPU time with large number of iterations, thus representing a challenging computational task.Table 5 compares the computational requirements by the two iterative schemes with five values of head tolerances.Since numerical oscillations occurred in the Newton scheme for 10 −2 m, this case will be excluded in the numerical discussion, and not shown in the plot (Figure 11) of convergence and MBEs.Simulations scenarios such as MBE, the total number of iterations used for each simulation, the average nonlinear iterations required for each step for the various simulations, the largest time stepping, and CPU time are included in this table.As shown in the table, the number of iterations for both iterative methods increases with the reduction of head tolerances.A perfect mass balance is guaranteed by the two schemes for all of the numerical experiments.We note here that 1147 and 46,165 steps are required to complete the simulation with head tolerance 10 −3 m for the Picard and Newton techniques, respectively.As compared with the Picard method, the Newton scheme required more than 46 times the iterations and 174 times the CPU time for simulating the 12-day drainage event.Furthermore, for the Newton method, the total steps and computational time increase severely with the reduction of tolerance.This is because discontinuity arises across the boundaries between different soil layers, and saturation gradients are steep on each side of the interface.The calculation of three derivative terms in the Jacobian makes the Newton scheme more costly; as a result, the numerical model needs to recompute (backstepping) many times than the current solution.A comparison of the behavior of convergence and MBEs for the Picard and Newton schemes is shown in Figures 11 and 12, respectively.The results are broadly representative, in that for most of the simulations, very little difference was found among the saturation predictions.Moreover, the MBEs usually follow parallel paths.Therefore, to avoid the excessive CPU time without significantly affecting convergence difficulties, one can use 10 −3 m of head tolerance.

Conclusions
A realistic computational approach is presented that effectively addresses the minimal head tolerance to solve the Richards equation accurately for four difficult test problems in combination with Picard and Newton linearization techniques.The tested examples and computational results presented illustrate clearly that the finite element formulation developed in this work is effective in handling severely nonlinear flow problems, including time-varying boundary conditions, nonhomogeneous soil hydraulic parameters, and an extremely sharp infiltration front in space.Simulations of redistribution, infiltration, and drainage problems in a one-dimensional soil column were used to evaluate the performance of the two schemes.An adaptive time-stepping scheme was incorporated easily in the Newton and Picard iteration techniques.The performance efficiency of the study is compared on the basis of MBE, CPU, the total number of time steps, and the average number of nonlinear iterations required by the algorithm to achieve the convergent solution within the specified tolerance.Accuracy and convergence profiles were generated to demonstrate the behavior of the two iterative techniques over a wide range of head tolerances.The main characteristic of the present formulation is to determine the least value of the head tolerances for the Picard and Newton iterative methods to obtain a stable and accurate solution with perfect mass balance.It has been shown that the schemes are extremely simple to implement, and such a procedure can greatly enhance the convergence properties of the iterative solution process with the head tolerance 10 −3 m.The numerical assessment has shown that the Picard scheme requires less CPU time per iteration than the Newton technique, because it needs to evaluate the Jacobian and produces a nonsymmetric coefficient matrix, which requires more computational effort.The procedure implemented in this work is very straightforward, and we believe its simplicity makes it an attractive approach for numerical simulations of variably saturated-unsaturated flow problems, and can easily be extended in more complex flow cases.

Figure 1 .
Figure 1.Computed pressure head profile of the Picard (left) and Newton (right) schemes at time 21,600 s for various tolerances.

Figure 2 .
Figure 2. Comparison of convergence and mass balance errors of the Picard scheme for various tolerances.

Figure 3 .
Figure 3.Comparison of convergence and mass balance errors of the Newton scheme for various tolerances.

Figure 4 .
Figure 4. Computed pressure head profile of the Picard (left) and Newton (right) schemes at time 17,280 s for various tolerances.

Figure 5 .
Figure 5.Comparison of convergence and mass balance errors of the Picard scheme for various tolerances.

Figure 6 .
Figure 6.Comparison of convergence and mass balance errors of the Newton scheme for various tolerances.
. A vertical discretization of 0.0125 m is used.Constant head boundary condition ψ(0, t) = 0.0 at the bottom of the domain, and a time-dependent boundary condition ψ(10, t) = −10 1.0 − 1.01e −t at the top of the domain with hydrostatic equilibrium initial conditions ψ(z, 0) = −z are applied.The time-varying boundary condition yields a difficult two-front problem.The van Genuchten model was fitted to the water retention curves, and the soil properties θ s = 0.301, θ r = 0.093, α = 5.47/m, n = 4.264.0,and K s = 5.040 m/days were used.

Figure 7 .
Figure 7. Computed pressure head profile of the Picard (left) and Newton (right) schemes at time 17,280 s for various tolerances.

Figure 8 .
Figure 8.Comparison of convergence and mass balance errors of the Picard scheme for various tolerances.

Figure 9 .
Figure 9.Comparison of convergence and mass balance errors of the Newton scheme for various tolerances.

Figure 10 .
Figure 10.Computed saturation profiles of the Picard (left) and Newton (right) schemes at time (approximately) 12 days for various tolerances.

Figure 11 .
Figure 11.Comparison of convergence and mass balance errors of the Picard scheme for various tolerances.

Figure 12 .
Figure 12.Comparison of convergence and mass balance errors of the Newton scheme for various tolerances.

Table 1 .
Computational performances of the Picard and Newton schemes for various tolerances.MBE: mass balance error; CPU: central processing unit.

Table 3 .
Computation performances of the Picard and Newton schemes for various tolerances.

Table 4 .
Soil hydraulic properties used in Test Problem 4.