Robustness and efficiency of iteration schemes for variably saturated flow across the range of soils, initial and boundary conditions found in practice

Water flow in partially saturated porous media is of importance in many disciplines such as hydrology, agriculture, environmental management or geotechnical engineering. A robust and accurate solution methodology applicable across the range of soils, initial and boundary conditions found in practice is difficult to identify. Three mass conservative iteration schemes for the approximation of the flow equation in partially saturated porous media are evaluated considering large contrasts in soil texture, initial moisture and boundary conditions typically found in real world application. The three schemes reviewed are the well established modified Picard algorithm (Celia-Scheme), and two promising approaches the L-Scheme and the Casulli-Scheme. A systematic variation of soil texture and initial moisture in an intense rainfall infiltration scenario is set up in order to consider the numerical challenges associated with high non-linearity. The present analysis provides insights into robustness and efficiency in solving variably saturated flow across the range of soils, initial and boundary conditions found in practice. Casulli’s scheme, in particular, provides an acceptable compromise between robustness and execution duration even under the severest conditions investigated, while the other schemes sacrifice one for the sake of the other. The authors believe that advances in the development of numerical schemes for partly saturated flow require the kind of methodology presented here in combination with relevant benchmark problems in order to enable a fair evaluation and comparison of numerical techniques across the related disciplines.


Introduction
Variably saturated flow problems appear in many disciplines such as hydrology, agriculture, environment and geotechnical engineering, e.g., in the assessment of agrochemical transport in soils, the estimation of surface run-off, the design of earthen dams and many other problems.The quasilinear partial differential equation attributed to Lorenzo A.
for variably saturated porous media.Characteristic issues encountered in real life applications are addressed considering heterogeneity, large contrasts in soil texture, varying initial moisture content and large intensity infiltration rates.All three schemes seek to enhance the robustness of the iterative solution compared to Newton's scheme.The well established modified Picard algorithm proposed by Celia et al. [3] represents the reference in this paper.Although it has been shown to suffer many of the same problems as Newton's scheme in engineering problems involving coarse soils [11,13], it rivals Newton's scheme in the frequency of employment in scientific literature and commercial software.The second scheme is the L-scheme originally proposed by Pop et al. [16,22].We follow the adaptation of the L-Scheme outlined in [11].L-Scheme was proven to be unconditionally convergent under the assumption of constant hydraulic conductivity (applicable only for nearly saturated or very fine grained soils) and is proven to be convergent even in the degenerative case (saturated or very dry conditions) [17].It should be noted that further developments for the L-Scheme have been proposed, for example the approaches by Mitra and Pop [14] or Stokkes et al. [19], which effectively switch locally between L-Scheme and either Celia's or Newton's method, respectively, depending on the current local state.We decided to focus on the original L-Scheme due to comparability reasons.The third scheme, proposed by Casulli and Zanolli [2], did not yet receive much consideration, at least in the engineering community.This scheme seems promising since its unconditional monotonicity and robustness for any time step size was mathematically proven and its performance for heterogeneous problems was demonstrated.Astonishingly, no comparison with other linearization schemes in terms of performance has been provided as of now.
There are various techniques with potential to improve the numerical solution of Richards' equation.For instance, Anderson acceleration can be applied to intermediate solutions, as demonstrated in [12] for Celia's scheme.Newton's method can be made more reliable through the use of globalization techniques such as line search and trust region method [8,10,20] or the combination with a more robust scheme [11].Another effective technique is the use of Kirchhoff transformation [7], which changes the primary variables to a more regular form of the problem.Similarly, numerical experience shows that different formulations of Richards' equation are favorable in terms of non-linearity in certain regions of soil moisture conditions, leading to the development of sophisticated variable switching techniques to harvest these advantages selectively [4,6].Additionally, nonlinearity tends to decrease with smaller time-step sizes, making adaptive time-stepping a widely used technique [9].In this paper, we focus on iteration schemes applied to the mixed form of the Richard's equation in order to provide a clearer understanding and direct comparison of these schemes.By avoiding adaptivity in time-stepping we intend to provide insights on the plain role of the iteration scheme.
In the following sections, we will present the governing equations for variably saturated flow and introduce the key aspects of each of the examined iteration schemes.In the third section, we will demonstrate the convergence properties of the considered schemes based on comparison with a onedimensional quasi-analytical solution.We will then illustrate a representative case considering rainfall infiltration, which has been identified to produce severe numerical problems under circumstances found in real world application.In the fifth section, we will compare the three iteration schemes in terms of efficiency and robustness.We show that Casulli's scheme, in particular, provides an acceptable compromise between robustness and execution duration for this kind of problem.Finally, we will provide an overview of our findings and give recommendations in relation to the investigated schemes in the sixth section.

Governing equations
In the framework of variably saturated flow, fluid mass conservation is expressed by the so-called mixed form of Richards' equation: In Eq. 1, ψ = 1/γ w ( p w − p a ) is the pressure head, with γ w the specific weight of water, p a and p w the pore-air pressure and pore-water pressure, respectively.θ(ψ) represents the volumetric water content, k u denotes the unsaturated hydraulic conductivity, which may vary with location x and reflects the dependence of hydraulic conductivity on water content.z is the height above a reference level (positive against the gravitational direction).The gas phase is assumed to be continuously connected to an atmospheric boundary and therefore always under atmospheric pressure ( p a ≡ 1 atm).Negative values of pressure head can be viewed as (negative) matric suction.Equation 1 demands soil hydraulic relationships.Water content θ(ψ) is assumed to be a nonlinear function of the pressure head (soil water characteristic curve or water retention curve).The unsaturated hydraulic conductivity is typically modeled as a function of either water content or pressure head.These soil hydraulic relationships are constitutive equations, of which the van Genuchten-Mualem model Eq. 2 is broadly applied.Constitutive models with fewer parameters like the Brooks & Corey model apply as well for the approach proposed by Casulli et al. [2] as will be shown later.
The model parameters must be estimated from direct measurements.θ r and θ s are the residual (ψ → − inf) and the saturated (ψ → 0) water content, respectively, α can be attributed to the air entry pressure, n and m = 1 − 1 n are shape parameters which describe the non-linear relationships between saturation, hydraulic conductivity and pressure head.k is the hydraulic conductivity under saturated conditions.Characteristic relationships are schematically displayed in Fig. 1.
In this study discretization of Eq. 1 is accomplished by means of cell-centered finite volumes in space using a two point flux approximation [21] and the Euler backward scheme in time.Regardless of the particular spatial discretization utilized, Eq. 1 is transformed into a system of nonlinear equations: Here, t denotes the time index and = ψ N i=1 is the vector of unknowns for the N finite volume cells.The Matrix Fig. 1 Schematic functional relationships of the van Genuchten -Mualem constitutive model H ( ) = θ i (ψ i ) I, with the N × N Identity Matrix I, represents the discrete water contents, F is the Flux Matrix containing the coefficients for the discretized Darcian fluxes due to pressure gradients and the right-hand side (RHS) comprises the discretized gravitational fluxes and known water contents from the previous time step.As mentioned above, the Newton method is generally applied to linearize the nonlinear terms in the matrices H and F. However, the locality of the Newton method strongly constraints the initial solution approximate which may render the approach impractical for certain conditions.Celia et al. [3] proposed an algorithm in which the hydraulic conductivity term is linearized by a Picard fixed-point iteration while the water content is linearized by a linear Taylor-Series expansion similar to Newton's method.Setting t,0 = t−1 as initial guess leads to the linearized system of equations with t and n denoting the time and iteration index respectively: The Jacobian C( ) = c i (ψ i ) I = ∂ ψ θ i (ψ i ) I may be interpreted as the ability of the soil to store or release water in the voids and is known as water capacity.Although Celia's modified Picard scheme is more robust than Newton's method, it still may fail under certain conditions or needs adaptive time-stepping to ensure convergence.Using the monotonicity of the θ(ψ), the L-Scheme transforms the iteration procedure to a contraction type [14]: The Jacobian C( ) is approximated by a fixed value L provided by the condition Eq. 5a.For the general case of variable hydraulic conductivity List and Radu [11] identified the optimal choice as L = L θ , which is used in the linearized system Eq.5b.
Casulli and Zanolli [2] identified the convex shape in the water capacity function c(ψ) (see Fig. 1) as a source for convergence problems of algorithm Eq. 4 and proposed a Jordan decomposition of c(ψ) into two monotonic increasing functions c(ψ) = p c (ψ) − q c (ψ) as shown in Fig. 2. Consequently, by independent integration (see Eq. 6) of the decomposed water capacity subfunctions, the water retention curve is decomposed into two functional relationships θ = θ 1 − θ 2 , as well (s.Fig. 2).
After decomposition of θ the system of equations becomes Keeping in mind that Q c represents the Jacobian for θ 2 , Newton's method applied only to the second term and, for the sake of readability, omitting the time and iteration indices (t, n) Eq. 7 is expressed with the outer Casulli iterates k := t,n,k as in Eq. 8a.
The initial guess 0 = n−1 is the solution of the previous Picard iteration.Linearizing the first term accordingly yields the nested Newton-type algorithm Eq. 8b with inner iterates k,l := t,n,k,l , with the solution of the previous outer iteration k,0 = k−1 as the initial guess.The nested iterations introduced in this scheme are terminated when condition Eq. 9a for the first and condition Eq. 9b for the second level of iteration are fulfilled.The right hand sides represent the L 2 norms of the inner and outer residuals.The tolerances ε p and ε q are each set to 10 −5 as proposed by Casulli and Zanolli [2].
This linearization has the advantage that p c and q c are monotonically increasing over the whole range.Interestingly, Casullis's scheme requires only monotonicity and not necessarily steadiness which makes it a good choice for simpler parametric functions as proposed by Brooks and Cory [1] or any characteristic curve composed of several branches with discontinuous slopes.

Choice of solver parameters
For all algorithms, convergence is accepted, if at the end of the outermost iteration step (n-loop) condition Eq. 10 is met.
Where, r is the L ∞ norm of the residuals of Eq. 3 with the current solution approximation ˜ t,n , φ being the vector containing the cell porosity and ε is the solution tolerance.ε can be interpreted as the mean saturation change accepted between two subsequent iterations.We examine suitable values for ε in the next section on basis of a one-dimensional homogeneous test case.We capped the maximum number of iterations per timestep at 1000.After this we let the solver continue with the next time step but mark the run time to signal convergence problems.All calculations have been carried out using the Open-FOAM v2012 library [21] on a AMD EPYC 7742 (Rome) (64 Core -2,25 Ghz) using one core per calculation.

Convergence properties
In this section, we consider a one-dimensional infiltration problem based on the quasi-analytical solution provided by [15].The test case involves a semi-infinite soil column with an initially uniform soil suction of ψ = −1m.The model domain extents over 10 m in depth, which ensures that the bottom boundary is sufficiently distant as not to affect the infiltration front.At the beginning of the calculation the soil suction at the top boundary is instantaneously increased to ψ = −0.05m(near full saturation).To assess the convergence properties of the three iteration schemes, we calculate the moisture distribution after 6h of infiltration and we compare how well they agree with each other and to the quasi-analytical solution.The parameters for the soil constitutive model (van Genuchten -Mualem) are listed in Table 1.The spatial domain is discretized using cells of 1cm length, and the time step is set to t = 100 s.The solution tolerance ε is varied in the range 10 −1 to 10 −5 (see Fig. 3).
In our study, we found that all of the solvers tended towards the same solution with decreasing values of ε.This solution did not exactly match the analytical solution at the infiltration front, due to numerical diffusion caused by spatial and temporal discretization.Yet, with only a couple of centimeter of deviation the accuracy was deemed as being sufficient.We observed that the different schemes did not all converge towards the analytical solution in the same way.While Celia's scheme tended to overestimate the infiltration depth and produced less stable solutions when the solution tolerance was set to the lowest setting, the L-Scheme underestimated the depth of the infiltration front under the same conditions.This follows from the constant term L, which is added to the main diagonal of the system matrix in the L-Scheme, which acts similarly to an increased storage or decreased time step, within the nonlinear iteration loop and is disappearing with convergence.The behavior of Celia's scheme can be similarly interpreted, while this scheme locally decreases or increases the storage, leading to overshoots as well as undershoots, respectively.For the smallest tolerance of ε = 10 −1 , Celia's and L-Scheme cannot be expected to converge, since they consist of only one set of iterations for which this tolerance is observably too high.In contrast, Casulli's scheme already converged to the settled solution with a tolerance of ε = 10 −1 , without noticeable difference when the tolerance was tightened to ε = 10 −3 .This is a consequence of the two nested iterations incorporated within this scheme.By the time the outer iteration is checked for convergence, the additional terms (acting as storage terms) introduced by the scheme have already vanished.These results suggest that the convergence of hydraulic conductivity is a peripheral concern as soon as water content is approximated with high precision, since this term is updated after the inner loops have converged.
Table 1 Hydraulic properties in the one-dimensional test case Fig. 3 Moisture fronts as water content θ over infiltration depth z after 6 hours

Real-world benchmark infiltration problem
To evaluate the strengths and weaknesses of the three different iteration schemes, a realistic 2D test-case from geotechnical engineering practice was set up.The problem considers infiltration during 72 hours due to precipitation (rainfall intensity i = 30 mm/d) into the backfill next to the foundation of a building.The analysis aims to predict the occurrence of saturated zones or the rise of groundwater level during or following the extreme rainfall event.These might require waterproof sealing and additional protection against uplift of the adjacent structure.The geometry and boundary conditions are presented in Fig. 4. The back-fill consists of coarse-grained sand placed in the excavation between the structure and the fine-grained in-situ soil.The mesh consists of unstructured triangular finite volume elements with an adequate resolution to represent the expected gradients (Fig. 4).The groundwater table (ψ = 0) is kept fixed at the bottom boundary.

Soil hydraulic properties
A series of calculations was performed for backfill materials of varying hydraulic properties covering the typical range of textures for this kind of rather coarse-grained material.In any soil hydraulic constitutive model the change of water content  2. A plot visualizing the parametric curves over pressure head for the examined backfill soils is presented in Fig. 5 for comparison.The value of α can be interpreted as a measure related to the number of coarse voids that can not hold water by capillary forces tightly enough to prevent gravity from pulling that water volume downward.The impact of α on the nonlinearity of the constitutive relationships is reflected in the drop of water content and the associated sharp decline in Table 2 Soil hydraulic properties of the backfill material according the van Genuchten -Mualem parametrization In  hydraulic conductivity over orders of magnitude within small changes of pressure head / water content.The assignment of the initial moisture content distribution represents an important issue in the analysis of infiltration dynamics.Furthermore, due to the contrast of the soil hydraulic properties, the initial water content of both backfill and in-situ soil was expected to affect the solver performance.Hence, the initial condition was set according to the steadystate solution of six different infiltration scenarios listed in Table 3.This choice considers initial moisture conditions in the range of very dry to very wet.

Size of timesteps
The choice of time step not only influences the accuracy but also impacts the solver performance.Solvers based on implicit time discretization must negotiate a balance between the size of time step that can be taken and the difficulty of the resulting system.To describe variably saturated flow a linearized system of equations must be solved, which tends to be easier to solve when using smaller time step sizes.As a result, it may be more efficient to solve multiple, well-conditioned systems with smaller time steps rather than fewer, larger time steps with ill-conditioned systems.The influence of nonlinearities also tends to be smaller over small time steps.
To study the influence of time step size all algorithms were tested with two different time step sizes, 1 h and 1/4 h.

Influence of soil hydraulic properties on the infiltration process
Infiltration is controlled by both soil hydraulic properties and antecedent water content distribution.6 illustrates the impact of these factors on infiltration dynamics in three different scenarios.In Fig. 6A, the calculation starts with no preceded precipitation (corresponding to a hydrostatic pressure distribution).In comparison, Fig. 6B corresponds to a precipitation of 400 mm/a, resulting in rather wet initial conditions and faster infiltration in both in-situ soil and backfill material.The specific discharge distributions (bottom row) in both cases show that the infiltration front is concentrated at the interface between backfill material and in-situ soil.This is due to the lower hydraulic conductivity of the sloped in-situ soil, which acts as a barrier and causes the moisture content at the interface to increase and hence the hydraulic conductivity to augment.In Fig. 6C, the backfill material has a substantially coarser texture than the previous cases, as indicated by the large contrast in the air-entry parameter α from 1.8 to 6.0 m −1 .The infiltration dynamics in the backfill material exhibit different characteristics, with a much shallower depth of the infiltration front.However, the infiltration front penetrates the coarser backfill material in a sharper fashion and leads to higher pressure heads, as seen in the pressure head variations with time at point P (Fig. 7).To summarize, a decrease in initial water content delays the infiltration front, and coarser soils (as indicated by larger air-entry parameters α) tend to delay the infiltration process but generate sharper infiltration fronts, leading to increased peak pressure.The study showed that, contrary to previous assumptions, the formation of a saturated zone adjacent to the building is unlikely, even for long duration rain events and high antecedent moisture levels.

Performance
In this section, we compare performance of the three schemes under consideration based on run-time costs using performance maps (see Fig. 8).These display CPU-time ratios τ on basis of the shortest execution times of any algorithm for a given test case.Based on that each value of τ in the performance plot shows how many times longer the respective scheme needed to finish a corresponding test case compared to the fastest algorithm for that problem.The value τ =1.0 means that an algorithm solved that particular test case the fastest.The underlying execution times correspond to 120 h of simulated time from start of the precipitation, since within this time period the maximum moisture increase from infiltration occurred (see Fig. 7).In Fig. 8 an asterisk behind the value of τ reveals that the iterative scheme did not converge to the required tolerance within 1000 iterations.Strictly from that point onward the solution might be inaccurate or even wrong.Inspection of these cases showed that the excess of the prescribed number iterations occurred only during a couple of time steps at the beginning of the computation when the gradients resulting from the prescribed rainfall intensity at the upper boundary are maximal.In the present study the solver was allowed to continue the computation as if convergence had been achieved.It is strictly stated that this is a risky strategy and it is up to the user to verify that the solution is physically based.Within this investigation with three different schemes such a verification is straight forward, however an evaluation may result impossible without any reference of what a plausible solution would be.
Based on the results, Celia's scheme showed a notable level of efficiency compared to the other two schemes, finishing fastest for almost all test cases it converged in.However, the scheme slightly loses efficiency when applied to coarser soils and drier conditions, which were also the conditions in which it showed weaknesses in terms of robustness.These kinds of problems are well known and were expected since imbibition into rather dry porous media generally implies strong gradients and therefore large variations of the storage and conductivity properties within small distances.Higher α values correspond to coarser soils representing higher degrees of non-linearity in the constitutive model which therefore affects the convergence rate of the scheme.Out of 30 test cases, the five cases with driest conditions and two cases with coarsest soils failed to converge, leading to a success rate of 80%.If full convergence is enforced (convergence within the allowed iterations for every time step), two more cases with α = 6 m −1 failed, diminishing the success rate to 73%.When the time step size was decreased to t = 1/4h, the scheme generally took longer to solve the test cases.The increased overhead for solving numerous time steps cannot be recovered and the algorithm turns out to be slower.With decreased time step size, only one more case was brought to completion, although only if full convergence is not enforced.A closer look reveals that even for some cases in which Celia's scheme with larger time step size converged fully, the smaller time step size hinders full convergence.This is a consequence of the harsh but realistic Fig. 6 Saturation distribution and isobars (top row), specific discharge rate q and isopotential lines (bottom row) for three sets of inital moisture conditions and soil textures at 42 hs after initiation of the rainfall event boundary condition, with instantaneously increased infiltration rate at the beginning of the calculation.Larger time step sizes allowed for this discontinuity to be "smeared out" by allowing more inflow to happen during the larger time interval prescribed by the time step size.A refined treatment of the boundary condition, such as ramping up the infiltration rate over time or refining the mesh density at the top boundary, are well-known strategies that might help mitigate this problem but introduce further complexity into the modeling task.
The L-Scheme was considerably slower compared to Celia's scheme, taking around 4 times as long even for the least critical conditions.Efficiency dropped further for drier initial conditions and coarser materials, where run-times ascended up to almost 20 times the run-time of the fastest algorithm.In terms of robustness, the L-Scheme displayed some advantage as all test cases could be brought to com-pletion by not enforcing full convergence.Yet, the success ratio without this loosened restriction was 66%, which was lower than for Celia's scheme.Decreasing the time step size for L-Scheme had a similar effect as for Celia's scheme, Fig. 7 Time variation of pressure head ψ at point P (see Fig. 6A).Lines are calculated with Casulli's scheme, circles with Celia's scheme and crosses using L-Scheme Fig. 8 Performance maps of CPU-time ratios τ over all test cases for each solver and time step size (left column: 1 h; right column: 1/4 h) with run-times generally increasing and previously fully converging test cases struggling.In fact, the impact of higher non-linearities associated with drier initial conditions and coarser materials was most pronounced for L-Scheme with a reduced time step size, as is predicted by the theory of the contraction rate of the linearization error given in [11].
Casulli's scheme was generally not the fastest even under moderate conditions (wetter initial conditions, less coarse).Yet, it was only between 20-30% slower compared to Celia's scheme.For wetter initial conditions, the 1 h-Casulli scheme almost equates in overall performance the 1/4 h-Celia scheme, while being always convergent.However, the reduced time steps of the Celia algorithm deliver the additional benefit of smaller local time-discretization truncation errors (although not necessarily global errors).Casulli's scheme displayed increased efficiency for both drier initial conditions and coarser materials, and the scheme took the lead for the toughest test cases.Casulli's scheme converged for every test case even without loosening the restrictions on allowed number of iterations, making it the most robust for this particular set of test cases.A closer look at the performance of the Casulli-Scheme with a 1/4 h time step reveals that Casulli's algorithm does not benefit from smaller time steps in this regard either.Unlike the other two schemes though, a decreased time step size did not increase the need to loosen the iteration restrictions for any test case.Casulli's scheme with a time step size of 1/4h converged for every test case.Nonetheless, the run-time increased (between 2-3 times), leaving increased accuracy associated with smaller time step sizes as the only advantage.

Conclusion
More than a few of recently published advanced schemes to the solution of Richards' equation were validated for rather moderate hydraulic soil properties and under somewhat academic initial and boundary conditions.In hydrology, agriculture or geotechnical engineering practice more often than not a much broader range of soil textures (typically quite coarse) under a variety of moisture conditions are encountered that cause high degrees of non-linearity when it comes to the solution of Richards' equation.Thus, to form a comprehensive basis to decide whether to adopt a promising scheme, comparative studies must take conditions into account that are relevant to the practice.
Based on two test cases we compared the well-established Celia's scheme [3] to two advanced iteration schemes, the L-Scheme [11,16] and Casulli's scheme [2].The choice of the solution tolerance is generally not independent of the particular iteration scheme.Aiming at a fair comparison of the three schemes a one-dimensional quasi-analytical solution was introduced to identify the respective solution tolerances for equivalent accuracy.Then, in a series of variably saturated flow computations of a representative problem comprising realistic conditions with respect to soil textures, initial soil moisture conditions and infiltration rates, the performance of all three schemes was investigated.Celia's Scheme displayed notable efficiency but showed the well known constraints for dryer conditions and coarser soils.The L-Scheme was considerbly slower compared to Celia's scheme, but managed to bring all test cases to completion.Casulli's scheme was slightly slower than Celia's scheme for moderate conditions, but displayed robustness and efficiency for all test cases including the driest initial moisture conditions and coarsest soil textures.The results show that for Celia's scheme as well as for the L-Scheme a decreased time step size did not enhance robustness significantly.
From the perspective of hydrological, environmental or geotechnical engineering practice, in which infiltration scenarios in rather coarse-textured soils and dry initial conditions are quite frequent the compromise between robustness and execution duration provided by Casulli's scheme is acceptable.Overall, the results provide insight into the trade-off between accuracy, efficiency and robustness when choosing an iteration scheme for variably saturated flow under the real world conditions.With the present analysis, the authors provide a benchmark test problem for variably saturated flow across the range of soils, initial and boundary conditions found in practice.The authors believe that advances in the development of numerical schemes for partly saturated flow require a repository of this kind of benchmark problems permitting a fair evaluation and comparison of numerical techniques across research groups.

Fig. 4
Fig.4 Considered geometry, spatial discretization and boundary conditions for the test cases (q n := surface normal flux).Backfill material in dark grey and in-situ soil in light grey

Fig. 5
Fig. 5 Water retention curves (top), water capacity (middle) and relative hydraulic conductivity (bottom) for the considered backfill soil textures corresponding to α-Parameters from 1.8 to 6 m −1 and up to a pressure head of ψ=-3 m

Table 3
Infiltration intensities q init for the steady state computations of the initial moisture distribution