A study on iterative methods for solving Richards` equation

This work concerns linearization methods for efficiently solving the Richards` equation,a degenerate elliptic-parabolic equation which models flow in saturated/unsaturated porous media.The discretization of Richards` equation is based on backward Euler in time and Galerkin finite el-ements in space. The most valuable linearization schemes for Richards` equation, i.e. the Newtonmethod, the Picard method, the Picard/Newton method and theLscheme are presented and theirperformance is comparatively studied. The convergence, the computational time and the conditionnumbers for the underlying linear systems are recorded. The convergence of theLscheme is theo-retically proved and the convergence of the other methods is discussed. A new scheme is proposed,theLscheme/Newton method which is more robust and quadratically convergent. The linearizationmethods are tested on illustrative numerical examples.


Introduction
There are plenty of societal relevant applications of multiphase flow in porous media, e.g. water and soil pollution, CO 2 storage, nuclear waste management, or enhanced oil recovery, to name a few. Mathematical modelling and numerical simulations are powerful, well-recognised tools for predicting flow in porous media and therefore for understanding and finally solving problems like the ones mentioned above. Nevertheless, mathematical models for multiphase flow in porous media involve coupled, nonlinear partial differential equations on huge, complex domains and with parameters which may vary on multiple order of magnitudes. Moreover, typical for the type of applications we mentioned are long term time evolutions, recommending the use of implicit schemes which allow large time steps. Due to these, the design and analysis of appropriate numerical schemes for multiphase flow in porous media is a very challenging task. Despite of intensive research in the last decades, there is still a strong need for robust numerical schemes for multiphase flow in porous media.
In this work we consider a particular case of two-phase flow: flow of water in soil, including the region near the surface where the pores are filled with water and air (unsaturated zone). By considering that the pressure of air remains constant, i.e. zero, water flow through saturated/unsaturated porous media is mathematically described by Richards' equation which has been proposed by L.A. Richards in 1930 (see e.g. [Bear and Bachmat(1991)]). In (1), Ψ denotes the pressure head, θ the water content, K stands for the hydraulic conductivity of the medium and z for the height against the gravitational direction. Based on experimental results, different curves have been proposed for describing the dependency between K, θ and Ψ (see e. g. [Bear and Bachmat(1991)]), yielding the nonlinear • The design of a new scheme based on the L−scheme and Newton's method, the L−scheme/Newton method which is robust and quadratically convergent.
• Provide the theoretical proof for the convergence of the L−scheme for Richards' equation (without using the Kirchhoff transformation) and discuss the convergence of modified Picard and Newton methods. The analysis furnishes new insights and helps towards a deeper understanding of the linearization schemes.
The present paper can be seen as a continuation of the works [Celia et al.(1990)] and [Lehmann and Ackerer(1998)], and it is written in a similar spirit. We added in the study the L−schemes (including a new scheme combining Newton's method with the L−scheme) and we focus now on 2D numerical results (the two mentioned papers based their conclusions on 1D simulations). We present illustrative examples, with realistic parameters so that the computations are relevant for practical applications.
The paper is structured as follows. In the next section the variational formulations (continuous and fully discrete) of Richards' equation are presented together with the considered linearization schemes. In Section 3 we discuss the theoretical convergence of the methods. The next section concerns the numerical results. The paper is ending by a concluding section.

Linearization methods for Richards' equation
Throughout this paper we use common notations in the functional analysis. Let Ω be a bounded domain in IR d , d = 1, 2 or 3, having a Lipschitz continuous boundary ∂Ω. We denote by L 2 (Ω) the space of real valued square integrable functions defined on Ω, and by H 1 (Ω) its subspace containing functions having also the first order derivatives in L 2 (Ω). Let H 1 0 (Ω) be the space of functions in H 1 (Ω) which vanish on ∂Ω, and H −1 (Ω) its dual. Further, we denote by ·, · the inner product on L 2 (Ω) or the duality product H 1 0 (Ω), H −1 (Ω), and by · the norm of L 2 (Ω). L f stays for the Lipschitz constant of a Lipschitz continuous function f (·).
We consider to solve the Richards' equation (1) on (0, T ] × Ω, with T denoting the final time and with homogeneous Dirichlet boundary conditions and an initial condition given by Ψ(0, x) = Ψ 0 (x) for all x ∈ Ω. We will use linear Galerkin finite elements for this study, but the linearization methods considered can be applied to any discretization method. We restrict the formulations and analysis to homogeneous Dirichlet boundary conditions just for the sake of simplicity, the extension to more general boundary conditions being straightforward (the numerical examples in Section 4 involve general boundary conditions). The continuous Galerkin formulation of (1) reads as: Find Ψ ∈ H 1 0 (Ω) such that there holds for all φ ∈ H 1 0 (Ω), with e z := ∇z. Results concerning the existence and uniqueness of solutions of (2) can be found in several papers, e.g. [Alt and Luckhaus(1983)].
For the discretization in time we let N ∈ N be strictly positive, and define the time step τ = T /N , as well as t n = nτ (n ∈ {1, 2, . . . , N }). Furthermore, T h is a regular decomposition of Ω ⊂ IR d into closed d-simplices; h stands for the mesh diameter. Here we assume Ω = ∪ T ∈T h T , hence Ω is assumed polygonal. The Galerkin finite element space is given by where P 1 (T ) denotes the space of linear polynomials on any simplex T . For details about this finite element space and the implementation we refer to standard books, like e.g. [Knabner(2003)]. By using the backward Euler method in time and the linear Galerkin finite elements defined above in space, the fully discrete variational formulation of (2) at time t n reads as: for all v h ∈ V h . At the first time step we take Ψ 0 h = P h Ψ 0 ∈ V h , with P h : H 1 0 (Ω) → V h being the standard projection. We assume in the next that the fully discrete schemes above have a unique solution and we refer to [Arbogast(1990), Arbogast et al.(1993), Ebmeyer(1998), Pop(2002)] for a proof.
At this point, dealing with the doubly nonlinear character of Richards' equation due to the relations K(θ) and θ(Ψ) is essential. We will briefly present in the following the main linearization methods used to solve the nonlinear problem (4): the Newton method, the modified Picard method (called simply Picard's method below, when does not exist a possibility of confusion) and the L−schemes.
We denote the discrete solution at time level n (which is now fixed) and iteration j ∈ N by Ψ n,j h henceforth. The iterations are always starting with the solution at the last time step, i.e. Ψ n,0 h = Ψ n−1 h . The Newton method to solve (4) reads as: holds for all v h ∈ V h . Newton's method is quadratically, but only locally convergent. As mentioned above, although Ψ n,0 h := Ψ n−1 h might be an appropriate choice, failure of Newton's method can occur (see [Radu et al.(2006)] and the numerical examples given below).
The modified Picard method was proposed by [Celia et al.(1990)] and reads as: holds for all v h ∈ V h . The modified Picard method was shown to perform much better than the classical Picard method [Celia et al.(1990), Lehmann and Ackerer(1998)]. The idea is to discretize the time nonlinearity quadratically, whereas the nonlinearity in K is linearly approximated. The method is therefore linearly convergent. The method still involves the computation of derivatives and in the degenerate case might also fail to converge (see the numerical examples in Section 4). The L−method was proposed for Richards' equation by [Yong and Pop(1996), Slodička (2002), ] and it is the only method which exploits the monotonicity of θ. The L−scheme to solve the nonlinear problem (4) reads as: holds for all v h ∈ V h . To ensure the convergence of the scheme, the constant L should satisfy L ≥ L θ (=: sup Ψ |θ (Ψ)|) (see Section 3 for details). The L−scheme is robust and linearly convergent. Furthermore, the scheme does not involve the computation of any derivatives. The key element in the new scheme is the addition of a stabilisation term L Ψ n,j h − Ψ n,j−1 h , v h , which together with the monotonicity of θ will ensure the convergence of the scheme.
Remark 2.1. It is to be seen that in the case of a constant K, the methods (5) and (6) coincide. Moreover, if L is replaced by the Jacobian matrix in (7) one obtains again the modified Picard scheme (6).
Any of the linearization methods presented above leads to a system of linear equations for Ψ n,j h (more precisely, the unknown will be the vector with the components of Ψ n,j h in a basis of V h ). The derivatives of the water content and the hydraulic conductivity in case of the modified Picard scheme and Newton's method can be computed analytically or by a perturbation approach as suggested by [Forsyth et al.(1995)] and occurring integrals are approximated by a quadrature formula.
For stopping the iterations, we adopt a general criterion for convergence given by with the Euclidean norm · and some constants ε a > 0 and ε r > 0. The tolerances ε a and ε r in criterion (8) are both taken as 10 −5 in all numerical simulations in this paper as proposed by [Lehmann and Ackerer(1998)]. We refer to [Huang et al.(1996)] for possible improvements of the stopping criterion. The Newton method is the only method out of the proposed three which is second order convergent. Nevertheless, it is not that robust as the other, linearly convergent, methods. In order to increase the robustness of Newton's method one can perform first a few (modified) Picard iterations, this being the combined Picard/Newton scheme proposed in [Lehmann and Ackerer(1998)]. The Picard/Newton method is shown to perform better than both the Newton and the modified Picard method [Lehmann and Ackerer(1998)]. We propose in this paper also a combination of the L−scheme with the Newton method, the L−scheme/Newton method. The mixed methods are based upon the idea to harness the robustness of the L−scheme or the modified Picard scheme initially and to switch to Newton's method e.g. if is satisfied for δ a , δ r > 0, similar to criterion (8). However, an appropriate choice of the parameters δ a , δ r is intricate and heavily dependent on the problem, for which reason a switch of the method after a fixed number of iterations may be an alternative. As shown in Section 4 this new method incorporating the L−scheme seems to perform best with respect to computing time and robustness.

Convergence results
In this section we will rigorously analyse the convergence of the L−scheme and discuss the convergence of Newton's and modified Picard's method. We denote by the error at iteration j. A scheme is convergent if e n,j → 0, when j → ∞.
The following assumptions on the coefficient functions and the discrete solution are defining the framework in which we can prove the convergence of the L−scheme.

(A2) The function K(·) is Lipschitz continuous and there exist two constants
We can now state the central theoretical result of this paper.
Theorem 3.1. Let n ∈ {1, . . . , N } be given and assume (A1) -(A3) hold true. If the constant L and the time step are chosen such that (16) is satisfied, the L−scheme (7) converges linearly, with a rate of convergence given by Proof. By subtracting (4) from (7) we obtain for any v h ∈ V h and any j ≥ 1 By testing the above with v h = e n,j and doing some algebraic manipulations one gets or, equivalently By using now the monotonicity of θ(·), its Lipschitz continuity (A1), the boundedness (from below) and Lipschitz continuity of K(·), i.e. (A2), the boundedness of ∇Ψ n h , and the Young and Cauchy-Schwarz inequalities one obtains from (13) After some obvious simplifications, the inequality (14) becomes Finally, by choosing L > 0 and the time step τ such that and by using the Poincare inequality (recall that e n,j ∈ H 1 0 (Ω)) e n,j ≤ C Ω ∇e n,j , from (15) follows the convergence of the scheme (7) We continue with some important remarks concerning the result above and the implications to the convergence of the Newton and modified Picard methods.
Remark 3.1. In the case of a constant hydraulic conductivity K (or if we refer to Richards' equation after Kirchhoff's transformation and without gravity, see e.g. [Pop et al.(2004)]), the condition for convergence of the L−scheme (7) simply becomes and there is no restriction on the time step size. Furthermore, the assumptions (A2) and (A3) are not necessary in this case.
Remark 3.2. The rate of convergence (11) depends on K m , τ and L, but it is independent of the mesh size. Smaller L or larger time steps are resulting in a faster convergence. We also emphasise that larger hydraulic conductivities will imply a faster convergence as well.
Remark 3.3. In the general case, the optimal choice is L = L θ and τ = , which is relatively mild because it does not involve the mesh size or any regularization parameter.
Remark 3.4. The convergence of the L−scheme is global, i.e. independent of the initial choice. Nevertheless, it is obviously beneficial if one starts the iterations with the solution of the last time step.
Remark 3.5. The convergence of the modified Picard method and of the Newton method is studied in [Radu et al.(2006)] for the case of constant K or for Richards' equation after Kirchhoff's transformation and without gravity. A regularization step is in this case necessary to ensure the convergence. The corresponding convergence condition to (16) will look like with denoting the regularization parameter and h the mesh size, d the spatial dimension and C a constant not depending on the discretization parameters. A similar condition is derived also for the Jäger-Kačur scheme, see [Radu et al.(2006)]. We remark that the condition (20) is much more restrictive than the condition (16). The proofs in [Radu et al.(2006)] are done for mixed finite element based discretizations, but the proof for Galerkin finite elements is similar. The condition (20) is derived by using some inverse estimates and it is in practice quite pessimistic. Nevertheless, we emphasise the fact that the convergence is ensured only when doing a regularization step (reflected by the in (20)) and this is what one sees in practice as well (see Section 4).
Remark 3.6. One can extend the convergence proof in [Radu et al.(2006)] for Newton's and modified Picard's methods to the general case of a nonlinear K and saturated/unsaturated flow. Under a similar assumption (A2) for the modified Picard and an assumption involving also the Lipschitz continuity of the derivative of K for Newton's method one can show the convergence of the methods. The modified Picard will be linearly convergent, whereas Newton's method will converge quadratically. The condition of convergence will be similar to (20) for both methods. From the theoretical point of view, only a quantitatively increased robustness for the Picard method comparing with Newton's method should be expected, i.e. when e.g. the mesh size becomes smaller if one of the method fails then also the other (see Fig. 2, where Newton's methods is not converging and modified Picard converges, but increasing the number of elements leads to divergence for the modified Picard or Picard/Newton methods as well). This is not the case with the L−scheme, which is clearly the most robust out of the considered methods, see Section 4.
Remark 3.7. By using error estimates derived as mentioned in the remark above, one can construct an indicator to predict the convergence of Newton's method. Based on this, one can design an adaptive algorithm for using the L−scheme only when necessary. Nevertheless, because the L−iterations are so cheap and the resulting linear systems are (much) better conditioned, it seems that the L−scheme/Newton is almost that fast as the Newton method. In Example 1 in Section 4 we even experienced that the L−scheme/Newton was faster than the Newton method. Therefore, we simply recommend the use of the L−scheme/Newton with a fixed number of L−iterations (4-5), without any indicator predictions. It the case of convergence failing, one should as a response automatically increase the number of L−iterations. We never experienced the need of more than 10 L−iterations in order to guarantee the convergence of the L−scheme/Newton.

Numerical results
In this section, numerical results in two spatial dimensions are presented. The considered linearization schemes: the Newton method, the modified Picard method, Picard/Newton, the L−scheme and the L−scheme/Newton are comparatively studied. We focus on convergence, computational time and the condition number of the underlying linear systems. We consider two main numerical examples, both based on realistic parameters. The first one was developed by us, the second is a benchmark problem from [Schneid(2004)]. Different conditions are created by varying the parameters. The sensitivity of the schemes w.r.t. the mesh size h is particularly studied. All computations were performed on a Schenker XMG notebook with an Intel Core i7-3630GM processor.

The relationships K(Ψ) and θ(Ψ) for both examples are provided by the van Genuchten-Mualem model, namely
in which θ S and K S denote the water content respectively the hydraulic conductivity when the porous medium is fully saturated, θ R is the residual water content and α and n are model parameters related to the soil properties. We compute the derivatives of K and θ analytically whenever they arise. The evaluation of integrals is executed by applying a quadrature formula accurate for polynomials up to a degree of 4.
Remark 4.1. The use of automatic differentiation might speed up the Newton method, but the concerns regarding the robustness will remain. This and the fact that most of the codes for solving Richards' equation do not have implemented automatic differentiation, were the reasons to compute the derivatives as mentioned above.

Example 1
This example deals with injection and extraction in the vadose zone Ω vad located above the groundwater zone Ω gw . The composite flow domain is Ω = Ω vad ∪ Ω gw defined as Ω vad = (0, 1) × (−3/4, 0) and Ω gw = (0, 1) × (−1, −3/4]. We choose the van Genuchten parameters α = 0.95, n = 2.9, θ S = 0.42, θ R = 0.026 and K S = 0.12 in parametrization (21). The choice n > 2 implies Lipschitz continuity of both θ and K. Constant Dirichlet conditions Ψ ≡ −3 on the surface Γ D = (0, 1)×{0} and no-flow Neumann conditions on Γ N = ∂Ω\Γ D are imposed. The initial pressure height distribution is discontinuous at the transition of the groundwater to the vadose zone and is given by Ψ 0 ≡ Ψ vad on Ω vad and Ψ 0 = Ψ 0 (z) = −z − 3/4 on Ω gw . We investigate two initial pressure heights in the vadose zone, Ψ vad ∈ {−3, −2}. In the vadose zone, we select a source term taking both positive and negative values given by f = f (x, z) = 0.006 cos(4/3πz) sin(2πx) on Ω vad , whereas we have f ≡ 0 in the saturated zone Ω gw . We examine the numerical solutions after the first time step for τ = 1. A regular mesh is employed, consisting of right-angled triangles whose legs are of length h = ∆x = ∆z for h ∈ { 1 10 , 1 20 , 1 30 , 1 40 , 1 50 , 1 60 } (the mesh size is actually h √ 2). The parameters regulating the switch for the mixed methods are taken as δ a = 2 and δ r = 0. The computation using the L−scheme was carried out with parameter L slightly greater than L θ = sup Ψ θ (Ψ) = 0.2341 for the given van Genuchten parametrization, to be specific L = 0.25. However, as pointed out in the analysis, when the influence of the nonlinear K is not that big (see Remark 3.1), a constant L bigger than L θ 2 is enough for the convergence. According to our experience, this is the limit relevant for the practice. Hence, we performed another computation with parameter L = 0.15. For the mixed L−scheme/Newton we choose L = 0.15 as well.
The results for Example 1 are presented in Figs. 1 -6 and discussed in detail below.

Convergence
In case of higher initial moisture in the vadose zone, that is Ψ vad = −2, convergence was observed for all methods and all investigated meshes. For the choice Ψ vad = −3, Newton's method failed on each mesh, the modified Picard scheme exhibited convergence only for h ≥ 1 40 , whereas both parametrizations of the L−scheme converged on all meshes. This is consistent with the theoretical findings in Section 3, in particular with Remark 3.6.

Numbers of iterations
The required numbers of iterations are depicted in Figs. 1 and 2. Missing markers indicate that the iteration has not converged. For either value of Ψ vad , the smaller parameter L = 0.15 in the L−scheme yielded the criterion for convergence to be fulfilled after fewer iterations than L = 0.25.
For Ψ vad = −2, the modified Picard scheme required less iterations than the L−scheme on coarse meshes, but for h ≤ 1 40 , it needed at least as many iterations as the L−scheme. Newton's method featured an even smaller number of iterations which was found to be independent of the mesh size in our computation. The number of iterations for the mixed Picard/Newton scheme did not differ significantly from the one for Newton's method, while the mixed L−scheme/Newton needed the least iterations.
For Ψ vad = −3, the modified Picard scheme had a benefit over the L−scheme in view of the number of iterations whenever it converged, although the number of iterations increased considerably as the mesh became finer. The mixed schemes gave the best results with respect to the number of iterations, the application of the mixed Picard/Newton scheme however being limited to coarse meshes.     Fig. 3 shows the simulation times for Ψ vad = −2. Although the modified Picard scheme needed less iterations than the L−scheme with L = 0.25, the differences of computation times were small, since the modified Picard scheme requires the computation of matrices including θ (Ψ). For Newton's method, K (Ψ) has to be calculated in addition. Nevertheless, it converged more rapidly than the modified Picard scheme. As reported by [Lehmann and Ackerer(1998)], combination of the modified Picard scheme and Newton's method further improved the performance in terms of computation time. However, both L−scheme with L = 0.15 and mixed L−scheme/Newton exhibited faster convergence than the mixed Picard scheme on dense grids, the mixed L−scheme/Newton only taking 71.2% of computation time compared to the mixed Picard/Newton scheme for h = 1 60 . The simulation times for Ψ vad = −3 are presented in Fig. 4. The mixed schemes computed the solution faster than the non-mixed schemes on each mesh, the mixed L−scheme/Newton taking roughly half the computation time in comparison to the non-mixed L−scheme with L = 0.15.

Condition numbers
In light of the accuracy of the numerical results, it is interesting to examine the condition numbers of the lefthand side matrices in the system of linear equations for the coefficient vector. Estimations for the condition numbers with respect to L 1 (Ω), denoted by · 1 calculated using the MATLAB function condest() are plotted in Figs. 5 and 6 for the non-mixed methods, averaged over all iterations. They did hardly differ from each other at several iteration steps and condition numbers for the mixed methods corresponded approximately to the ones of the respective non-mixed method in each iteration. For both values of Ψ vad , the L−scheme with L = 0.25 featured the lowest condition numbers, followed by its counterpart with L = 0.15. In case of Newton's method being convergent, it exhibited higher condition numbers than the L−scheme. In all computations, the condition numbers in the modified Picard scheme were the highest, furthermore, they increased most rapidly when it came to denser meshes. All methods required more iterations and computation time when the vadose zone was taken to be dryer initially and the arising matrices were worse-conditioned than for the moister setting.

Example 2 (Benchmark problem)
In order to compare the linearization methods in the numerical simulation of a recognized benchmark problem, we consider an example used by [Haverkamp et al.(1977)], [Knabner(1987)] and [Schneid(2004)] amongst others. It describes the recharge of a groundwater reservoir from a drainage trench in two spatial dimensions (Fig. 7). The domain Ω ⊂ R 2 represents a vertical section of the subsurface. On the right hand side of Ω, the groundwater The initial and boundary conditions are taken as in which n denotes the outward pointing normal vector. Initially, a hydrostatic equilibrium is thus assumed. The computations are undertaken for two sets of parameters adopted from [van Genuchten(1980)], characterising silt loam respectively Beit Netofa clay. For both soil types, the solution is computed over N = 9 time levels. The time unit is 1 day and dimensions are given in meters. The van Genuchten parameters employed as well as the parameter ∆t D governing the time evolution of the upper Dirichlet boundary, the time step τ and the simulation end time T are listed in Table 1. We used a regular mesh consisting of 651 nodes. The simulations invoking the L−scheme were carried out with L = sup Ψ θ (Ψ) (referred to as L−scheme 1) and with L slightly smaller (referred to as L−scheme 2) for both soil types, that is L = 4.501 · 10 −2 and L = 3.500 · 10 −2 for the silt loam soil and L = 7.4546 · 10 −3 and L = 6.500 · 10 −3 for the clay soil. The mixed methods switched to Newton's method when condition (9) held true for δ a = 0.2 and δ r = 0. All the considered linearization methods converged for both soil types. The pressure profiles computed with mixed L−scheme 2/Newton at time T are presented in Fig. 8 and are as expected for this benchmark problem. Table 2 shows the total numbers of iterations, the computation times and the average of the estimated condition numbers of the left-hand side matrices with respect to · 1 , in case of mixed methods split up in the two involved schemes. In what follows, the foregoing numerical indicators, i.e. the number of iterations, the computational time and the condition numbers are to be discussed in detail.

Numbers of iterations
As to the non-mixed methods, it is not surprising that more complex methods yielded smaller numbers of iterations, i.e. Newton's method converged after the fewest iterations, followed by the modified Picard scheme. L−scheme 2 had the edge over L−scheme 1, but still needed some more iterations than the modified Picard scheme for both soil types. The numbers of iterations of the mixed methods exhibit a salient result: The advantage of the modified Picard scheme over the L−scheme with regard to the number of iterations vanished when coupling the schemes to Newton's method and the mixed L−scheme 2/Newton required less iterations than the mixed Picard/Newton scheme. This suggests that the L−scheme stands out due to a rapid approach towards the solution in the first iteration steps. Among all methods, Newton's method provided convergence after the least number of iterations for both van Genuchten parametrizations.

Computation times
When it comes to the comparison of computation times, it is striking that the performances of the methods substantially varied between the simulations for silt loam and Beit Netofa clay. While Newton's method featured the shortest computation time among the non-mixed methods in case of silt loam owing to the low number of required iterations, computation in case of the clayey soil took long using Newton's method as compared to the L−scheme.
In the silt loam simulation, computation times of the L−scheme were clearly greater than the ones of Newton's method, but switching to Newton's method vastly improved the computation time so that the L−scheme 2/Newton turned out to be the fastest method. In contrast, the computation times for the clay soil demonstrate that in some cases, switching to Newton's method may even be disadvantageous. Although the mixed L−scheme/Newton converged in fewer iteration steps than the non-mixed ones, changing to Newton's method provoked a deterioration of the computation time. This might indicate that the L−scheme be less susceptible to parametrizations of the hydraulic relationships lacking of regularity than the modified Picard scheme and Newton's method since the hydraulic conductivity for the parametrization of the Beit Netofa clay is not Lipschitz continuous. The modified Picard scheme was found to be the slowest method for the silt loam soil, the computation time for Beit Netofa clay was barely less than the one related to Newton's method.

Condition numbers
In view of the condition numbers of the left-hand side matrices, the L−scheme excelled for both soil types: The condition numbers with either value of L were remarkably lower than the ones arising when Newton's method or the modified Picard scheme were employed, to be more specific by a factor of minimum 11 for the silty soil and still by a factor of minimum 5 for the clayey soil. Apparently, incorporation of the derivative of the water content entailed a considerable deterioration of the condition. The virtual equality of the condition numbers for the modified Picard scheme and Newton's method was probably due to the proximity of the solution to a hydrostatic equilibrium which caused the only term distinguishing Newton's method from the modified Picard scheme in equation (5) to be small because of ∇Ψ n h ≈ −e z .

Conclusions
In this paper we considered linearization methods for the Richards' equation. The methods were comparatively studied w.r.t. convergence, computational time and condition number of the resulting linear systems. The analysis was done in connection with Galerkin finite elements, but the schemes can be applied to any other discretization method as well, and similar results are expected. We focused on the Newton method, the modified Picard method, the Picard/Newton and the L−scheme. We proposed also a new mixed scheme, the L−scheme/Newton which seems to perform best. We conducted a theoretical analysis for the L−scheme for Richards' equation, showing that it is robust and linearly convergent. We also discussed the convergence of the modified Picard and Newton methods.
The L−scheme is very easy to be implemented, does not involve the computation of any derivatives and the resulting linear systems are much better conditioned as the modified Picard or Newton methods. Although it is only linearly convergent, seems to be not much slower than the Newton (or Picard/Newton) method, and in some cases even faster. The L−scheme is the only robust one, a result which can be shown theoretically and it is supported by the numerical findings. Only a relatively mild constraint on the time step length is required.  Furthermore, when the hydraulic conductivity K is a constant, there is no restriction in the time step size. In this case the only condition necessary for the global convergence of the L−method is L ≥ L θ 2 . We proposed a new mixed scheme, the L−scheme/Newton which is more robust than Newton but still quadratically convergent. This new mixed method performed best from all the considered methods with respect to computational time. Even in cases when Newton converges, the L−scheme/Newton seems to be worth, being faster for the examples considered.
The present study is based on two illustrative numerical examples, with realistic parameters. The examples are two dimensional. One of the examples is a known benchmark problem. The numerical findings are sustaining the theoretical analysis.