Ross scheme, Newton–Raphson iterative methods and time-stepping strategies for solving the mixed form of Richards’ equation

. The solution of the mathematical model for ﬂow in variably saturated porous media described by the Richards equation (RE) is subject to heavy numerical difﬁculties due to its highly nonlinear properties and remains very challenging. Two different algorithms are used in this work to solve the mixed form of RE: the traditional iterative algorithm and a time-adaptive algorithm consisting of changing the time-step magnitude within the iteration procedure while the nonlinear parameters are computed with the state variable at the previous time. The Ross method is an example of this type of scheme, and we show that it is equivalent to the Newton– Raphson method with a time-adaptive algorithm. Both algorithms are coupled to different time-stepping strategies: the standard heuristic approach based on the number of iterations and two strategies based on the time truncation error or on the change in water saturation. Three different test cases are used to evaluate the efﬁciency of these algorithms.Thenumerical results highlight the necessity of imple-menting an estimate of the time truncation errors.

Abstract. The solution of the mathematical model for flow in variably saturated porous media described by the Richards equation (RE) is subject to heavy numerical difficulties due to its highly nonlinear properties and remains very challenging. Two different algorithms are used in this work to solve the mixed form of RE: the traditional iterative algorithm and a time-adaptive algorithm consisting of changing the timestep magnitude within the iteration procedure while the nonlinear parameters are computed with the state variable at the previous time. The Ross method is an example of this type of scheme, and we show that it is equivalent to the Newton-Raphson method with a time-adaptive algorithm.
Both algorithms are coupled to different time-stepping strategies: the standard heuristic approach based on the number of iterations and two strategies based on the time truncation error or on the change in water saturation. Three different test cases are used to evaluate the efficiency of these algorithms.
The numerical results highlight the necessity of implementing an estimate of the time truncation errors.

Introduction
Water movement in soils is one of the key processes in the water cycle since it contributes to the renewal of groundwater resources through recharge, to vegetation growth through transpiration, to soil fertility through salinization/alteration and to atmospheric humidity through evaporation and transpiration. Water movement is usually modeled using the Richards equation (Richards, 1931), which is now commonly adopted for many studies in soil science and/or hydrology, including the use of physically based hydrological models applied to large-scale catchments and for long time simulations (e.g., for climate change studies). However, this equation is highly nonlinear, and despite numerous efforts over the last 40 years, its numerical solution requires much computational time.
Assuming a rigid solid matrix, the Richards equation (RE) is given by where θ is the volumetric water content (L 3 L −3 ), S w is the water saturation (-), s 0 accounts for fluid compressibility (L −1 ), ψ is the pressure head (L), q is the water flux based on the extended Darcy's law (L T −1 ), t is the time (T ), z is the vertical coordinate (positive upward) (L), f is the sink/source term (T −1 ), K is the saturated hydraulic conductivity tensor (L T −1 ) and k r (ψ) is the relative hydraulic conductivity (−). The model includes initial and boundary conditions of the Dirichlet (prescribed pressure head) or Neumann (prescribed flux) type.
Equation (1) is also called the mixed form of RE. Two alternative formulations of the mixed form exist for RE.
The pressure form is defined by Published by Copernicus Publications on behalf of the European Geosciences Union.
Constitutive relations are required to solve RE. For the pressure-water content relationship, the most common model is the Van Genuchten model (van Genuchten, 1980): where m = 1 − 1/η, S w is the effective saturation, θ r and θ s are the residual and saturated volumetric water content, respectively, and α and η are experimentally estimated coefficients. This model is usually associated with the Mualem model (Mualem, 1976) for the relative permeability of the aqueous phase: A summary of the most popular relations can be found in Belfort et al. (2013).
Adaptive time-stepping strategies based on time truncation error control were found to be superior to other approaches (Hirthe and Graf, 2012;Kavetski et al., 2001;Tocci et al., 1997). The method of lines using the DASPK integrator was applied to the Richards' equation by Matthews et al. (2004), Miller et al. (1998), andTocci et al. (1997), among others. The method of lines consists of discretization of the spatial part of the PDE only, leading to a system of ordinary differential equations. It has been found to be significantly more efficient than other temporal discretizations (Miller et al., 2006). However, Kavetski and Binning (2002b) reported difficulties in obtaining convergence for the DASPK solver associated with an arithmetic mean of inter-block conductivities for the most difficult problem addressed by Miller et al. (1998). Additionally, very few non-iterative schemes have been developed Binning, 2004, 2002a;Paniconi et al., 1991).
Despite the many existing numerical methods, solution of RE is still a challenging research topic, with many remaining questions about reduction of the computational time, treatment of nonlinearities, and improvement of the accuracy of these methods for difficult problems such as infiltration in very dry soils (Diersch and Perrochet, 1999;Forsyth et al., 1995;Hills, 1989). The need for efficient algorithms to solve this equation has increased during recent decades because it has been recognized that explicit modeling of flow in the unsaturated zone has to be implemented in land surface models (Vergnes et al., 2012). In their recent review of land surface models, Clarke et al. (2015) push for a mechanistic modeling of the flow in soils. They consider the implementation of the mixed form of the Richards equation to be an improvement to the modeling of soil moisture variations. They also underline the need for efficient algorithms to solve RE, to allow the implementation of stochastic approaches and/or automatic parameter estimations.
In this study, we analyzed the performance of different algorithms based on the Newton-Raphson method since the classical Picard scheme has been found to be less efficient (Lehmann and Ackerer, 1998). Applied to the soil moisture form of RE, we demonstrate that the recently developed Ross method (Ross, 2003;Crevoisier et al., 2009;Zha et al., 2013) is equivalent to the Newton-Raphson method (Sect. 2). A detailed presentation of the Newton-Raphson method applied to the mixed form or RE is given in Sect. 3. The standard Newton-Raphson algorithm is based on the computation of the corresponding matrices in an iterative way by updating the parameters until convergence. An alternative algorithm has been suggested more recently where the parameters are kept unchanged within one time step and the time step is adapted to reach convergence. This algorithm has been applied to the pressure-based form of RE by Kavetski and Binning (2002a) and to the soil moisture form by Crevoisier et al. (2009), Ross (2003, and Zha et al. (2013). Although this algorithm is called "non-iterative" because the parameters are not updated during the calculation, iterations may be necessary to adapt the magnitude of the time step. Therefore, in the following, we will refer to the usual algorithm as "iterative" and to the alternative algorithm as "time-adaptive". To our knowledge, this alternative algorithm has never been applied to the mixed form of RE. Section 4 is dedicated to both algorithms and to the time-stepping strategy used for solving RE. Finally, in Sect. 5, the numerical accuracy and robust-

The Ross method and the Newton-Raphson method
The moisture-based formulation is applicable in unsaturated conditions only and is prone to numerical difficulties in the case of heterogeneous soils, explaining the reduced attention directed to this formulation. However, discontinuous water content can be handled by adapted schemes, and the moisture-based formulation appears to be very accurate for initially dry conditions (Zha et al., 2013(Zha et al., , 2015. Ross (2003) suggested a non-iterative formulation that has been recently extended to different soil conditions (Crevoisier et al., 2009;Varado et al., 2006a) and to two and three dimensions (Zha et al., 2013).
In its initial one-dimensional finite-volume formulation and for a volume (cell) i, the Ross method (Ross, 2003) is based on the following set of equations: where S n+1 i is the water saturation at cell/node i at time (n+1), q σ − (or q σ + ) is the water flux between cell i and (i − 1) (or i+1) at time t = t n + σ t, σ ∈ [0, 1] and z is the size of cell i. θ s,i is the saturated water content and θ r,i is the residual water content. For simplicity, we assume here that all cells are of the same size.
The previous mass balance Eq. (6) leads to the following equation for cell i: The Newton-Raphson method was initially developed as a root-finding algorithm of an arbitrary equation that has been generalized for solving a system of nonlinear equations. Applied to the soil moisture form of RE and using an implicit scheme, the NR consists in defining a residual based on the mass balance equation (Eq. 6) at iteration k for time step n+1 and for cell i written as where R n+1,k i is called the residual. The NR consists in computing the solution at iteration k+1 by estimating the residual of the next iteration R n+1,k+1 i using a first-order Taylor development and setting it equal to zero as The derivatives of this residual are which leads to the following set of linear equations: For the first iteration, we have S n+1,k+1 Whatever the formulation of the fluxes q (as a function of the pressure (see Eq. A1) or the water content, expressed by Kirchhoff transform as in Ross, 2003, or not), the implicit Ross method (Eq. 8 with σ = 1) is equivalent to the first iteration of the Newton-Raphson method (Eq. 13).

Newton-Raphson method for the mixed form Richards' equation
Because the pressure-based formulation does not ensure mass conservation -except for the approximation provided by Rathfelder and Abriola (1994) -and due to the limitations of the moisture-based formulation (see previous section), the mixed formulation has been widely used since the work of Celia et al. (1990). The mixed form of the Richards equation given by Eq.
(1) is rewritten as and is discretized by where A is the discretized form of the divergence term, B and E are the discretized forms of the storage terms, F is the discretized form of the sink/source term and the boundary conditions, n is the time step and k is the iteration counter. t n+1 is the time-step magnitude defined by t n+1 = t n+1 − t n . Matrices A, B, and E and vector F depend on the numerical scheme used for the spatial discretization. The implicit scheme is applied for the spatial discretization.

Algorithms and time-stepping strategy
The usual algorithm used to solve RE consists in defining a time step that remains constant and in iteratively computing the parameters and variables in the following way.
For a given time step n Define the time-step length t n+1 depending on the time-stepping strategy.
2. Computation of the system matrix R and the residual R.
4. Check convergence. If convergence is achieved, exit.
enddo Next time step where k is the iteration counter and maxit the maximum number of iterations.
The time-adaptive algorithm consists in calculating the nonlinear parameters with the pressure heads computed at time step n and adapting the time-step length. The algorithm is described by the following.
For a given time step n Computation of the variable θ n , the parameter K n and their derivatives dθ n dψ n , ∂K n ∂ψ n using ψ n . do k = 1, maxit 1. Define a time step t n+1,k depending on the timestepping strategy.
2. Computation of the system matrix R and the residual R.
4. Check convergence. If convergence is achieved, exit.
enddo Next time step The main advantage of the alternative algorithm is its avoidance of the computation of the variable θ , the parameter K and their derivatives dθ dψ and ∂K ∂ψ during the iterations. Due to the highly nonlinear relations between θ , K, dθ dψ , ∂K ∂ψ and the pressure, this computation may require significant CPU time.
The most popular time-step management during the simulation is that of the heuristic type (Miller et al., 2006). The Hydrol. Earth Syst. Sci., 21, 2667-2683, 2017 www.hydrol-earth-syst-sci.net/21/2667/2017/ time step t n+1 is computed depending on t n and the number of iterations k necessary to reach convergence in the following way: where k 1 , k 2 , m 1 , and m 2 are user-defined constants.
Other heuristic time-step management procedures have been suggested by Kirkland et al. (1992) based on the water volumes exchanged between the adjacent cells of the grid, and by Ross (2003), where the time-step size is controlled by the maximum allowed change in the saturation.
For the Ross method, the fluxes are computed first and the time-step magnitude is calculated accordingly using where S max is the user-defined maximum allowed saturation change. After the computation of the change in the saturation S, the time step is modified if the maximum of the computed change exceeds and the system of equations is solved again. More details about handling the fluxes at boundaries and saturated conditions can be found in Crevoisier et al. (2009), Ross (2003 and Varado et al. (2006b). The adaptive scheme used in this work evaluates the time steps through truncation error due to the temporal discretization as proposed by Thomas and Gladwell (1988). This scheme was already applied to the pressure-based formulation by Kavetski et al. (2001) and to the moisture-based formulation by Kavetski and Binning (2004).
The difference between the first-order and second-order time approximations can be considered as an estimate of the local truncation error of the first-order scheme. The firstorder approximation is given by The second-order approximation is using ∂ψ n+1 ∂t = ∂ψ n ∂t + t n+1 ∂ 2 ψ n ∂t 2 .
This truncation error is given by when the truncation error is smaller than γ , the temporal truncation error tolerance defined by the user, and the size of the next time step calculated by When the truncation error is larger than γ , the computation is repeated with a reduced time step defined as follows: where r max and r min are user-defined constants used to avoid overly drastic changes in the time step. s is considered to be a safety factor that ensures that the time-step changes are reasonable. EPS is used to avoid floating point errors when the truncation error becomes too small.

Evaluation of the algorithms' performance
We applied the NR method to the mixed form of RE using the standard iterative algorithm and the time-adaptive algorithm. A cell centered finite-volume scheme for the spatial discretization with an implicit Euler scheme for the temporal discretization has been used to solve the partial differential equation and arithmetic means are used to compute the interblock hydraulic conductivity. The detailed discretizations of the matrix R (ψ n+1,k ) and the vector R(ψ n+1,k ) (see Eq. 18) are given in Appendix A. The time-adaptive algorithms have been applied as described by the authors: Ross (2003) for the time stepping based on the saturation changes and Kavetski et al. (2001) for the time stepping based on the truncation errors.
For the standard iterative algorithm, we defined two types of errors to check the convergence: the error based on the maximum change in the state variables between two iterations defined by ε ψ = max i ψ n+1,k+1 i − ψ n+1,k i and the truncation error ε t defined by Eq. (24). Convergence is assumed to be achieved when where τ ψ,a and τ ψ,r are the absolute and relative user-defined tolerances and ψ n+1,k+1 imax is the pressure corresponding to ε ψ  and when where τ t,a and τ t,r have the same meaning as those for the previous criterion but ψ n+1,k+1 imax represents the pressure value corresponding to ε t .
The tested algorithms are summarized in Table 1. Computations of all possible combinations for the standard iterative scheme have been performed. We present only the four most efficient algorithms. We also analyzed convergence based on the nonlinear residual. It was found to be less restrictive than the previous criteria. Due to the definition of the NR method, the residual tends to zero, but it does not ensure a small value of ε ψ . Therefore, the results related to the reduction of the nonlinear residuals are not reported.
We investigated three one-dimensional problems with various initial and boundary conditions and hydraulic functions to assess the accuracy, efficiency and computational costs of the different algorithms. The selected test cases represent a range of difficult infiltration problems widely analyzed in the literature.
-TC1: infiltration in a homogeneous initially dry soil with constant prescribed pressure at the surface and prescribed pressure at the bottom (Celia et al., 1990); -TC2: infiltration in a homogeneous soil initially at hydrostatic equilibrium with a prescribed constant flux at the soil surface and prescribed pressure at the bottom (Miller et al., 1998); and -TC3: infiltration/evaporation in an initially dry heterogeneous soil, with variable positive and negative fluxes at the surface and free drainage at the base of the soil column (Lehmann and Ackerer, 1998).
For the three test cases, the soil hydraulic functions were described by Mualem-Van Genuchten models (Mualem, 1976;van Genuchten, 1980); see Eqs. (4) and (5). The required parameters, boundary conditions and initial conditions are summarized in Table 2. The evolution of the relative hydraulic conductivity, the water saturation and the specific moisture capacity with respect to the pressure values are shown in Figs. 1, 2 and 3, respectively. For TC1, the pressure will vary from −1000 to −75 cm only due to the specific conditions of this test case. Therefore, the parameter variations are smaller than those for the other test cases. Since the parameters' variations are more abrupt for test cases 2 and 3, their solutions are more challenging.
Preliminary tests were performed to define the optimal spatial discretization; i.e., a finer spatial discretization provided very similar results for a given convergence criterion and a given time-stepping strategy. Therefore, we can assume that the errors only originate from the time-step size and the linearization.
The following criteria were used for the time-stepping strategy: k 1 = 0.80, k 2 = 1.20, m 1 = 5, and m 2 = 10, which are the usual values for the heuristic strategy defined by Eq. (19); and Hydrol. Earth Syst. Sci., 21, 2667-2683, 2017 www.hydrol-earth-syst-sci.net/21/2667/2017/ Table 2. Domain size (L), initial conditions (IC), boundary conditions at the soil surface (BC u ) and at the soil bottom (BC l ), saturated hydraulic conductivity (K s ), residual and saturated water contents (θ r , θ s ) and shape parameters (α, η) for the different test cases. K M (t) is the hydraulic conductivity of the last grid cell. Length and time units are centimeters and seconds, respectively. r min = 0.10, r max = 4.0, s = 0.9, and EPS = 10 −10 , which are the standard values for the time-stepping scheme based on the time discretization error defined by Eq. (26) (Kavetski et al., 2001).
To perform a consistent comparison of the time-stepping strategies, the maximum allowed change in saturation (see Eqs. 20 and 21) has been evaluated using the maximum change in the pressure, according to the following relationship: The simulations have been performed using different values of τ r and with τ a = 0.0. We used several criteria to evaluate the performance of these codes. A typical error used in solving RE is the global cumulative mass balance error defined by where z i is the size of the cell/element i, θ n+1 i is its water content at time t n+1 , θ 0 i is the initial water content, and q k in and q k out are the inflow and outflow, respectively, at the domain boundaries at time t k . M is the number of cells/elements. The fluxes at the boundaries are defined by q k = 1 2 q k + q k−1 . The mass balance errors were checked for each run but were found to be negligible since we solved the mass-conserving RE form.
While it is necessary to satisfy the global mass balance for an accurate numerical scheme, a low mass balance error is not sufficient to ensure the accuracy of the solution. Therefore, solutions have also been compared with the reference solution obtained using a very fine temporal discretization and the iterative Newton-Raphson method. This comparison is based on the average relative error defined by where M is the number of cells, ψ ref is the reference solution andψ is the tested numerical solution. ε 1 represents the average absolute relative error (called L 1 -norm in the following), ε 2 is the average quadratic error (L 2 -norm) and ε ∞ is the highest local relative difference between the numerical and reference solutions (L ∝ -norm).
Since the time-adaptive algorithm does not require the computation of the parameters and their derivatives during the iterative procedure, we use N sol to denote the number of times where the system of equations is solved and N param to denote the number of times where the parameters are computed. Of course, these counters are equal to each other for the standard algorithm, which leads to computational costs depending on 2N sol . N param is less than N sol for the timeadaptive algorithm. For comparison purposes, the computational costs are estimated by N sol for the standard algorithm and by (N sol +N param )/2 for the time-adaptive algorithm. The efficiency of the algorithms has been evaluated by comparing the computational costs for a given relative tolerance τ r . The errors are presented in the tables and the figures. The figures show some additional results not listed in the tables that already contain much information.

TC1: Infiltration in a homogenous soil with constant boundary conditions
This test case simulates an infiltration into a homogeneous porous medium. This problem is addressed here because it has been widely analyzed previously by many authors like Bouchemella et al. (2015), Celia et al. (1990), El Kadi and Ling (1993), Rathfelder and Abriola (1994), and Tocci et al. (1997), among others. The computations were performed with a spatial discretization of 0.1 cm. The initial time-step size was set to 1.0 × 10 −5 s, and the maximum time-step size was set to 400 s. The results for the iterative and time-adaptive algorithms are presented in Tables 3 and 4, respectively. When both convergence criteria are used (algorithms SH_ ψ_ t and SS_ ψ_ t), N trunc represents the number of times where the truncation error is the most restrictive condition. For the heuristic time-stepping schemes, the convergence is mostly linked to the truncation error (N trunc is close to N sol ), whereas when the saturation time-stepping scheme is used, the most restrictive criterion is the maximum difference in the pressure.
When the time-stepping scheme is based on saturation, for both iterative and time-adaptive algorithms, the number of iterations required to solve the problem is proportional to the relative tolerance. Therefore, highly accurate solutions incur high computational costs.
For the time-adaptive scheme, the number of parameter changes N param is close to the number of iterations for low tolerance values. Small tolerance values lead to small time steps, avoiding time-step adjustments. This is not the case for larger tolerance values that lead to larger time steps and therefore to additional iterations (see for example TA_T for the tolerance of τ r = 10 −2 - Table 4).
The three types of errors provide the same information. The best solution for one type of error is also the best solution for the other two.
On average, the iterative algorithm is faster than the timeadaptive algorithm that requires more iterations for a given error. This is also shown in Fig. 4 that presents the convergence rate of the L 2 -norm with respect to the computational costs, i.e., the number of iterations or number of iterations and number of parameter changes. The time-adaptive algorithm with time stepping based on the truncation errors performs quite poorly compared to the other algorithms. Irrespective of the tolerance, this algorithm leads to a wetting front moving faster (Fig. 5).
When the relative tolerance is set to a very low value (τ r = 10 −5 ), the iterative scheme with time stepping based on the saturation changes shows behavior that is different from that found for the less restrictive tolerance. The criterion based on truncation errors is no longer significant (N trunc = 252), possibly explaining why the accuracy of the scheme remains constant. This also indicates that errors due to time discretization have to be handled, either in the convergence criterion or in the time-stepping strategy.
For this test case, the most efficient algorithms are the iterative algorithms using the time-stepping strategy based on truncation error (ST_ ψ) or based on the saturation changes (SS_ ψ_ t). Saturation-based time-stepping strategies (SS_ ψ_ t and TA_S) show a linear decrease Hydrol. Earth Syst. Sci., 21, 2667-2683, 2017 www.hydrol-earth-syst-sci.net/21/2667/2017/ in L 2 with computational costs. For very high precision (L 2 < 10 −4 ), ST_ ψ outperforms the other algorithms. No convincing explanation has been found for the insignificant change in accuracy for SS_ ψ_ t at high precision.

TC2: Infiltration in a homogenous soil with hydrostatic initial conditions
This test case models an infiltration in a 200 cm vertical column of unconsolidated clay loam with non-uniform grain size distribution and was considered by Miller et al. (1998) to be a very challenging test. This problem was found to be more challenging from the numerical point of view compared to TC1 due to the relative permeability function that enhances the nonlinear behavior of Richards' equation (Figs. 1,  2, and 3). The cell size has been set to 0.125 cm, the initial time step to 10 −5 s and the maximum time-step magnitude to 1000 s. The different norms for the iterative and time-adaptive schemes are given in Tables 5 and 6.
Investigation of this test case leads to similar qualitative conclusions when the time-stepping scheme is based on the saturation differences (SS_ ψ_ t and TA_S). The standard scheme SH_ ψ fails to provide an accurate solution within a reasonable number of iterations (less than 10 7 ).
The most efficient methods are the schemes using the timestepping strategy based on truncation errors (Fig. 6). However, as found for TC1, the adaptive time algorithm TA_T failed to provide highly accurate results (L 2 -norm error less than approximately 4.5 × 10 −4 ). Figure 7 shows the time-step magnitudes for approximately equal L 2 -norms for the two time-adaptive algorithms and for the iterative algorithm using truncation errors www.hydrol-earth-syst-sci.net/21/2667/2017/ Hydrol. Earth Syst. Sci., 21, 2667-2683, 2017 Table 5. Relative errors and number of iterations obtained for the iterative algorithm depending on different convergence criteria for TC2 (n.c.: non convergence in less than 10 7 iterations).
Algorithm for time stepping (4.254 × 10 −4 within 3503 iterations for ST_ ψ, 4.563 × 10 −4 within 3098 iterations for TA_T and 4.844 × 10 −4 within 11358 iterations for TA_S). The timestep evolution is very similar for the three strategies: a linear increase until around 0.1 s, followed by a very slow increase until 20-30 s and a regular increase until the end of the simulation. ST_ ψ and TA_T strategies lead to the same time steps when the time reaches 1 s. The time-step sizes remain smaller for TA_S, which explains the significantly higher number of iterations required to solve this test case.

TC3: Infiltration/evaporation in a heterogeneous soil
This case study simulates infiltration in an initially dry heterogeneous soil with a succession of rainfall and evapora-tions as upper boundary conditions during 35 days. This problem differs from the two previous cases by the soil heterogeneity and also by the non-monotonic boundary conditions at the soil surface. It is expected that non-monotonic discontinuous boundary conditions will increase the difficulty in finding accurate solutions. The soil profile consists of three 60 cm thick layers. The layers are discretized using cells with the size of 0.10 cm. The prescribed fluxes are changing every day. For a given time, these fluxes are linearly interpolated. To avoid an overly rough time discretization of these boundary conditions, the maximum time-step magnitude has been fixed at 0.20 days. The initial time step is set to 10 −5 days. The relative errors estimated by the iterative algorithms and the time-adaptive algorithms are presented in Tables 7  and 8, respectively, and are plotted in Fig. 8.
Hydrol. Earth Syst. Sci., 21, 2667-2683, 2017 www.hydrol-earth-syst-sci.net/21/2667/2017/ Table 7. Relative errors and number of iterations obtained for the iterative algorithm depending on different convergence criteria for TC3 (n.c.: non convergence in less than 10 7 iterations, * convergence failed for 10 −3 , τ r = 0.90 × 10 −3 ). The standard iterative scheme fails to converge within the maximum number of iterations (10 7 ) when the tolerance is not sufficiently restrictive. The detailed analyses of the computation showed that the time-step size was quite large compared to the more restrictive conditions until day 28.0, where the infiltration fluxes were equal to 1.50 cm day −1 and where the conditions were near saturation due to the previous infiltration period. This led to a decrease in the time step to close to the minimum value (10 −8 s), causing the procedure to stop. More restrictive conditions lead to smaller time steps from the beginning of the simulation and a better approximation of the solutions during the entire simulation.
The iterative scheme coupled with the truncation-based time-stepping strategy showed surprisingly unstable behavior for τ r = 10 −3 . The scheme did not converge for τ r ∈ 0.96 × 10 −3 ; 1.04 × 10 −3 . The results presented in Ta-ble 7 and Fig. 8 are obtained for τ r = 0.90 × 10 −3 . At this stage of our work, we were not able to provide a meaningful explanation for this effect. The time-adaptive algorithm with the saturation-based time-stepping scheme is the most efficient for an L 2 -norm greater than 10 −4 . For more accurate results, the iterative method with the time-stepping strategy using the truncation error must be preferred. The impact of the time-stepping strategy for these two algorithms is shown in Fig. 9 for approximately the same L 2 -norm (2.051 × 10 −3 within 1283 iterations for TA_S and 1.517 × 10 −3 within 6504 iterations for ST_ ψ). The time-step changes are related to the boundary condition variations, as expected. The strategy based on the saturation variation leads to a longer time step than the strategy using the time truncation error. This difference can be quite important (see the simulation between days 25  and 30). The consequences of this difference are a reduced number of iterations but also a less accurate computation, irrespective of the error norm.

Summary and conclusions
The solution of RE is complex and very time-consuming due to its highly nonlinear properties. Several algorithms have been tested for the mixed form of the Richards equation, including time-adaptive methods. Based on the numerical examples that differ in their parameters (level of nonlinearity) and in their initial and boundary conditions, the conclusions and recommendations are the following.
1. Our numerical developments showed that the method suggested by Ross (2003) in its implicit formulation can be considered as a Newton-Raphson method with a time-adaptive algorithm.  2. The different algorithms have different convergence rates (accuracy improvement of the scheme as a function of the computational costs). Therefore, an algorithm can be very efficient for a given accuracy and less efficient for another level of precision. However, for these three test cases and, on average, the best performance in terms of efficiency was obtained using a stopping criterion based on truncation error with its corresponding time-step strategy (ST_ ψ). Similar results were obtained by Kavetski et al. (2001) for the pressurebased RE and by Kavetski and Binning (2004) for the moisture-based RE.
3. The mass balance is not a good criterion for the evaluation of the results because the mixed form preserves the mass balance, irrespective of the pressure distribution within the profile.
Hydrol. Earth Syst. Sci., 21, 2667-2683, 2017 www.hydrol-earth-syst-sci.net/21/2667/2017/ 4. The time truncation error should be implemented in numerical codes using the standard iterative procedure. The use of the maximum variable difference between two successive iterations only, which is usually implemented, does not provide any information about the accuracy of the time derivative approximation.
Our one-dimensional examples showed that the timeadaptive algorithm TA_T is very sensitive to the type of problem to be solved. The time-adaptive algorithm TA_S was less efficient than the usual schemes. However, for a larger number of elements like in two-dimensional or three-dimensional problems, this conclusion might be different because the time dedicated to the computation of the parameters can be significantly higher, unless tabulated values are used to evaluate the parameters and the required derivatives. Depending on the type of the problem that must be solved (parameter behavior with respect to the pressure, time variations of the boundary conditions), the time truncation errors may be predominant compared to the error corresponding to the pressure changes between two successive iterations. Therefore, we recommend the implementation of this stopping criteria associated with the time-stepping strategy as defined by Kavetski et al. (2001).

Appendix A
The numerical method used in the paper is implicit standard finite difference. For a cell i of the grid, the unsaturated flow Eq. (4) can be discretized in the following way: where n is the time step, K i− is the inter-block conductivity between cell i and (i − 1) defined by K i− = z i−1 K(ψ i−1 )+ z i K(ψ i ) z i−1 + z i , and K i+ is the inter-block conductivity between cell i and (i+1) defined by K i+ = . z i− = 1 2 ( z i−1 + z i ) is the distance between the center of cell (i − 1) and i. z i+ = 1 2 ( z i + z i+1 ) is the distance between the center of cell i and (i+1).
The residual is where k is the iteration counter.
In the case of prescribed flux at the upper boundary, the residual is written as Using the derivatives as defined in Eqs. (A5) and (A6), the matrix coefficients are changed as follows: If the pressure is described at the top of the soil, the corresponding flux is defined by and the derivative is The corresponding residual and the matrix coefficients are R 1 = z 1 θ n+1,k 1 − θ n 1 + S w s 0 ψ n+1,k 1 − ψ n 1 + t q n+1,k 1+ − q n+1,k The numerical code is written in FORTRAN 90 and is available upon request.