1 Introduction
Hydrogeological modeling is widely used for problems such as water resources management and safety assessment of hazardous objects and requires numerical modeling of groundwater flow. The objects of interest are often located near the surface, in partially saturated zone. In these cases, groundwater flow is described by Richards equation richards1931capillary ; bear2010modeling . Nonlinearity of the Richards equation limits applications of analytical methods and creates significant difficulties in numerical solution process, which are mainly attributed to arising nonlinear systems. Efficient solution techniques for both transient and steadystate cases are subject of extensive research farthing2003efficient ; farthing2017numerical .
This paper focuses on improvement of a special solution approach for boundary value problems for the steadystate Richards equation, namely, the nonlinearity continuation method anuprienko2021nonlinearity . In this procedure, the governing Richards equation in original problem is parametrized in a way that its ”degree of nonlinearity” can be controlled. Varying the parameter, one can get a sequence of problems from a simple linear problem to the original one. Each problem in the sequence is then solved with an iterative solver (Newton, Picard or their combinations anuprienko2021comparison ) with the solution of previous problem given as initial guess. Following allgower2003introduction , this procedure belongs to the class of predictorcorrector continuation methods with a trivial, zerothorder predictor step and the corrector step being application of an iterative nonlinear solver.
This paper contributes to the nonlinearity continuation method for the steadystate Richards equation in the following ways.

The nonlinearity continuation method is formulated in predictorcorrector terms, the firstorder predictor is applied and its effect on performance is examined;

Discretization framework, which was previously limited to finite volume schemes, is extended to mimetic finite difference method in order to test the nonlinearity continuation procedure on structurally different nonlinear systems.
The paper is organized as follows. The second section contains mathematical description of the groundwater flow model, the third section describes discretization and resulting nonlinear systems, the fourth section formulates nonlinearity continuation in predictorcorrector terms, the fifth section presents results of numerical experiments, and the conclusion provides some discussion.
2 Groundwater flow problem in variably saturated conditions
2.1 Governing partial differential equations
From the mass balance equation for water
(1) 
and the Darcy law for water flux
(2) 
the steadystate Richards equation richards1931capillary ; bear2010modeling
(3) 
is derived. Here the following variables are used:

– water flux;

– capillary pressure head;

– volumetric water content in medium;

– hydraulic conductivity tensor, a 3
3 s.p.d. matrix; for isotropic medium this tensor is scalar and can be characterized by one number ; 
– relative permeability for water in medium;

– specific sink and source terms.
It is perhaps more convenient to work with the steadystate Richards equation (3) in terms of hydraulic head (here is the vertical coordinate). Indeed, the relative permeability can be rewritten as a function of hydraulic head in the following way: and for simplicity denoted as . Then, one can formulate steadystate Richards equation in the following form:
(4) 
The same applies to the flux expression (2), and it becomes
(5) 
To completely define the boundary value problem, boundary conditions and additional constitutive relationships between water content , water pressure head and relative permeability are required.
2.2 Boundary conditions
The following boundary conditions are used for the steadystate Richards equation:

Dirichlet: specified hydraulic head ;

Neumann: specified normal flux .
2.3 Van Genuchten – Mualem constitutive relationships
This model is widely used in vadose zone applications and is capable of predicting capillary effects in soil. It is based on nonlinear functions for and proposed by van Genuchten van1980closed and Mualem mualem1976new :
(6) 
(7) 
where is the water content at full saturation, is the residual water content, is the effective saturation and , and are parameters of the model.
3 Discretization techniques
Complex geometries and layered structure of geological domains, as well as presence of small objects like wells, imply use of unstructured grids in groundwater modeling applications. Two types of grids are considered in this paper: triangular prismatic grids with occasional occurrence of other polyhedra and hexahedral grids with octree refinement and possibility to cut cells for better representation of the boundaries plenkin2015adaptive . In practical applications, these meshes often contain flattened cells that may be some arbitrary polyhedra.
Therefore, the discretization scheme should be able to deal with cells of quite general shape. Other desirable properties include mass conservation on discrete level, absence of nonphysical solutions and low computational complexity. Discretizations considered in this paper include finite volume (FV) and mimetic finite difference (MFD) methods. These are techniques for discretization of diffusiontype operators and have been employed in various applications such as simulation of groundwater flow or multiphase flow in porous media. While finite volume schemes have been already employed in earlier work anuprienko2021nonlinearity ; anuprienko2021comparison , the mimetic scheme is considered for the first time within the nonlinearity continuation method and is included in order to test the solution strategy for discretization other than finite volume.
3.1 Finite volume schemes
The finite volume method (see, for example, eymard2000finite ) is a common choice for subsurface applications due to its local conservativity and ability to handle meshes with polyhedral cells. After integrating equation (4) over a cell and using the divergence theorem, one gets a sum of normal fluxes over the cell faces. Approximation of these fluxes is the key point of a finite volume scheme. The simplest linear twopoint flux approximation (TPFA) remains the most popular choice for its simplicity. However, TPFA does not give approximation for nonorthogonal grids, which are very common in practical applications. Other FV schemes include linear multipoint flux approximations aavatsmark1998discretization ; agelas2008mpfa ; agelas2010g and schemes with nonlinear two and multipoint flux approximations that are monotone or satisfy discrete maximum principle for diffusion equation le2005schema ; kapyrin2007family ; danilov2009monotone ; misiats2013second ; nikitin2016nonlinear ; terekhov2017cell .
Of all finite volume schemes described, the following schemes are considered in this work:

linear TPFA;

MPFAO scheme aavatsmark1998discretization ;

nonlinear twopoint flux approximation (NTPFAB) terekhov2017cell ;

nonlinear multipoint flux approximation (NMPFAB) terekhov2017cell .
The relative permeability for a face can be either calculated as a halfsum of values from the two neighboring cells (central approximation) or taken from the cell with greater water head value (upwind approximation). In the author’s experience, central approximation leads to better convergence of nonlinear solvers since the structure of arising system does not change during iterations, but results in saturation oscillations. Upwind approximation gives more physically reasonable saturation distributions, but can result in severe convergence problems of nonlinear solvers.
Discretization leads to a system of nonlinear equations which takes the form
(8) 
where is a matrix with solutiondependent coefficients.
3.2 Mimetic finite difference scheme
The mimetic finite difference method is a way to construct discretization schemes with discrete analogues of differential operators that mimic properties of their continuous counterparts. An extensive overview of MFD history may be found in lipnikov2014mimetic , while theoretical foundations, derivation of mimetic schemes and their applications are presented in da2014mimetic . Recently, mimetic schemes have been applied in subsurface flows simulation, including groundwater flow described by Richards equation lipnikov2016new ; coon2020coupling and multiphase flows abushaikha2020fully .
Mimetic approach benefits from Richards equation considered as two firstorder PDEs (1) and (5). These equations represent diffusion problem in mixed formwith nonlinear solutiondependent coefficient. The scheme used in this work chooses the divergence as the primary operator and uses finitevolumelike approximation of it in cells, while the approximation of the gradient operator is built in such a way that it is negatively adjoint to the discrete divergence in appropriate discrete spaces with their scalar products. Details on the discrete gradient construction can be found in lipnikov2016discretization . The presence of solutiondependent relative permeability in (5) makes the construction of the scheme more complicated, since some approximations of nonlinear coefficient may lead to incorrect solutions for transient equations lipnikov2016mimetic . Following lipnikov2016new , upwind approximation of relative permeability was used in present work.
The resulting mimetic scheme uses cellcentered hydraulic head unknowns (one value per cell) and facecentered flux unknowns (one value per face, representing the normal flux). The respective nonlinear system, therefore, has significantly larger size compared to cellcentered finite volume schemes. That complexity is the price to pay for a number of attractive properties, among which are

ability to work with a wide range of cell types, including degenerate and nonconvex cells;

simultaneous calculation of fluxes;

easy incorporation of boundary conditions.
Another important property of the scheme is its similarity to finite volume schemes, which results from the same assignation of mesh unknowns and the same approximation of divergence operator. This ensures local conservativity and allows for easy incorporation into existing finite volume framework, including some constitutive relationships and procedures designed specifically for finite volume paradigm anuprienko2018modeling .
4 Nonlinearity continuation method in predictorcorrector terms and first order predictor
4.1 General ideas
Both finite volume and mimetic schemes result in a system of nonlinear equations of the form
(9) 
where
is the respective vector of discrete unknowns (for finite volume schemes it consists of cellbased head unknowns, while for the mimetic scheme it also contains facebased fluxes). An obvious way to solve the system is to apply standard iterative solvers such as Newton and Picard methods. However, these fail often. Newton method needs a good initial guess, finding which for a complicated flow structure is a hard task which becomes harder in case of an advanced discretization. Picard needs contraction properties of discrete operator and often suffers from slow convergence or lack thereof. Additional difficulties result from upwinding of the relative permeability, since it can change function
during iterations. Relaxation and line search methods paniconi1994comparison ; farthing2003efficient ; armijo1966minimization , as well as using Picard iterations to improve the initial guess paniconi1994comparison ; anuprienko2021comparison , may improve the robustness, but still are not effective enough.A more powerful approach is the nonlinearity continuation method anuprienko2021nonlinearity . The relative permeability function , which makes the equation nonlinear, is parametrized with a continuation parameter such that the resulting function, , has the following properties:

;

.
From these properties it follows that with one gets the original steadystate Richards equation, and with one gets a linear equation representing the full saturation case.
Several options to choose function exist, for example:

linear function ;

power function .
In this paper, only power fucntion is considered.
Nonlinear systems arising from discretization of the problems for parametrized Richards equation with parameter are denoted as
(10) 
where the subscript in emphasizes that solution corresponds to a certain value of the continuation parameter.
The nonlinearity continuation method in its previously studied form anuprienko2021nonlinearity ; anuprienko2021comparison is a stepping procedure. It starts with and gradually increases it to . This way a sequence of problems with different is obtained. Each problem in the sequence is solved with an iterative solver: Newton, Picard or a mix of the two, all combined with a line search technique. The solver takes solution of the previous problem as initial guess. For the first problem, with , there is no such initial guess, and a constant hydraulic head distribution is chosen as initial guess; however, since this problem with is linear (it represents the fully saturated case), iterative solvers converge easier (for linear discretization schemes they converge in 1 iteration). The number of problems in the sequence or, in other words, number of steps from to , is not known a priori. These steps are chosen within the solution process.
4.2 Predictorcorrector formulation
The nonlinearity continuation procedure described above lies in the class of predictorcorrector methods allgower2003introduction , where the corrector step is application of an iterative solver and the predictor step is setting initial guess for the said solver. Outline of the method with emphasized predictor and corrector steps is presented in algorithm 1. In its form studied in anuprienko2021nonlinearity ; anuprienko2021comparison , the nonlinearity continuation method has a trivial, or zerothorder predictor: it simply sets solution of the previous problem as initial guess for the next problem.
4.3 Firstorder predictor
In algorithm 1 the predictor step is setting initial guess for . Previously, a trivial, zerothorder predictor was used, which has the form
(11) 
A more sophisticated predictor may be used. In this socalled firstorder predictor a special vector is calculated, which describes sensitivity of the solution to the continuation parameter . After a problem with is solved, this vector is calculated from the linear system
(12) 
where is the Jacobian matrix calculated for the last obtained solution, and the righthand side is approximated with a finitedifference expression. The value of of is often chosen to be around square root of machine epsilon knoll2004jacobian for rather simple problems, and in this work is chosen.
4.4 Details on the corrector
Different nonlinear solvers may be used as correctors. Previous work anuprienko2021comparison studied Newton and Picard solvers as well as their combination. Experiments showed that the use of mixed solver may help to complete the nonlinearity continuation procedure in 1 step for simpler constitutive model from anuprienko2018modeling , but requires setting parameters, such as initial number of Picard iterations and initial relaxation coefficient. For the VGM model, the mixed solver gave no clear advantage. This, and the fact that the firstorder predictor improves initial guess (which is more crucial for the Newton method) are the reasons to use pure Newton method in this paper. Additionally, line search in accordance with Armijo rule is used armijo1966minimization ; anuprienko2021comparison , but it is not applied for the first 5 iterations. Convergence is controlled by the tolerances and , and not more than iterations are allowed.
5 Numerical experiments
The purpose of numerical experiments is to study effect of firstorder predictor on the number of continuation steps and nonlinear iterations, compared to the zerothorder predictor.
5.1 Implementation details
All methods described in this paper are implemented in C++ as a part of the GeRa (Geomigration of Radionuclides) software kapyrin2015integral ; gerasite . GeRa was originally designed for safety assessment of radioactive waste repositories and nowadays is a generalpurpose hydrogeological modelling tool. Computational core of GeRa is based on INMOST vassilevski2020parallel ; inmostsite , a platform which provides tools for numerical modeling on general unstructured grids with builtin MPI parallelization.
Of the discretization schemes used, only TPFA and MPFAO are available in GeRa. NTPFAB and NMPFAB schemes are included for research purposes from an external package terekhov2017cell . The MFD scheme has been implemented by the author recently and has not been used extensively yet.
Automatic differentiation module of INMOST is used to obtain the Jacobian matrix, and the corresponding linear systems are solved with Inner_MPTILUC, which is an internal INMOST solver based on BiCGSTAB preconditioned by second order CroutILU decomposition konshin2015parallel
with inversebased condition estimation
li2003crout and maximum product transversal reordering duff1999design .5.2 Capillary barrier test
The first problem, the capillary barrier test, describes infiltration in a twolayered domain. It was first presented in oldenburg1993numerical , and later used for verification of modeling software, including TOUGH2 webb1997generalization , FEFLOW diersch2013feflow and GeRa kapyrin2017verification . The problem is essentially twodimensional and features domain (presented at figure 1) which is approximately 100 m long and 6 m high. The domain is composed of two inclined layers (angled at 5%), each 0.5 m thick: the upper layer is sand, the lower layer is gravel. The media parameters are given in the table 1. On the lower and right boundaries zero hydraulic head condition is imposed, on the top boundary the infiltration rate of 0.0048 m/day is prescribed.
An important feature of the capillary barrier is the high contrast in media properties. In fully saturated conditions gravel is much more permeable than sand. However, in partially saturated state it is less permeable, which leads to occurrence of the capillary barrier effect. The complicated structure of the total conductivity curve contributed to the choice of power continuation function in the present work.
, m/day  , 1/m  

Sand  18.144  3.9  5.74  0.154  0.39 
Gravel  8640  490  2.19  0.011  0.42 
Numerical experiments were conducted on a 3200cell hexahedral grid and the following Newton parameters: , . Upwind approximation of was used. Computed hydraulic head and water saturation distrbutions (which vary very slightly for different discretizations) are depicted in figure 2, while solution characteristics are presented in table 2. For the first three schemes the firstorder predictor was able to decrease total computation time by up to 40% by reducing number of continuation steps and Newton iterations. For the NMPFAB scheme, the firstorder predictor decreased the number of steps and iterations too, but the total time did not decrease since the number of function evaluations in the line search increased. For the MFD scheme, the firstorder predictor gave no advantage. It should be noted that solution with the MFD scheme was the fastest, although it is the most expensive scheme on periteration basis, and the reason for this is the lowest number of iterations.
, s  # of successful (failed) steps  # of Newton iterations  

TPFA, 0  89.7  33 (35)  829 
1  51.9  27 (28)  463 
MPFAO, 0  375.0  31 (33)  667 
1  281.5  24 (25)  432 
NTPFAB, 0  48.7  35 (37)  409 
1  29.3  26 (27)  265 
NMPFAB, 0  63.0  32 (34)  493 
1  63.2  25 (27)  413 
MFD, 0  35.8  13 (14)  117 
1  57.0  23 (24)  198 
5.3 Realistic case
The second problem involves groundwater flow in a vertical crosssection of a reallife waste repository located near ground surface. The domain (see fig. 3) is approximately 90 m long and 12 m high and is composed of 6 isotropic materials with their properties listed in table 3. Rainfall recharge rate of 0.0001 m/day is prescribed on the top boundary, while on other boundaries hydraulic head vale of 0.5 m is imposed.
Several features contribute to difficulty of solving this problem:

Hydraulic conductivity ranges from 0.048 to 100 m/day with jumps up to 3 orders of magnitude at some interfaces;

Relative permeability curves (7) are nonsmooth for some materials;

Dry zones as well as multiple fully saturated ones are present;

Central approximation of results in significant unphysical oscillations of the water saturation, and upwind approximation should be used instead, which leads to convergence problems;

As often done in practice, computational mesh was coarsened in larger subdomains to reduce number of cells. This resulted in considerable nonorthogonality and problems with TPFA were expected.
The mesh consisted of 2800 cells. Solution was performed with the following Newton parameters: , . The problem turned out to be troublesome with only two schemes finishing: TPFA and MFD. Other schemes failed due to inability to find admissible . Solution process characteristics are presented in table 4 and show that the firstorder predictor reduced total computational time by more than 50% for TPFA and by almost 40% for MFD. In MFD case, the total number of continuation steps increased while the total number of iterations decreased. This may be attributed to the fact that inner Newton solver converges faster when firstorder prediction is good and diverges faster when it is not.
Calculated water saturation distributions are presented in figure 4. As expected, TPFA produced some numerical artifacts, namely, ”checkerboarding” in saturation (see figure 5) which likely results from inconsistent flux approximation. MFD, on the other hand, produced reasonable results.
, m/day  

Clay  0.048  0.8  1.09  0.068  0.38 
Soil  0.1  5.9  3  0.1  0.35 
Gravel  100  6  3  0.04  0.3 
Clay loam  0.2496  3.6  1.56  0.078  0.43 
Sand  7.128  14.5  2.68  0.045  0.43 
Repository  0.1  14.5  3  0.045  0.4 
, s  # of successful (failed) steps  # of Newton iterations  

TPFA, 0  47.2  17 (36)  383 
1  20.1  11 (23)  186 
MFD, 0  335.5  18 (39)  510 
1  185.5  19 (41)  312 
6 Conclusion
The nonlinearity continuation method, which has been recently applied for boundary value problems for the steadystate Richards equation, can be considered as a predictorcorrector procedure. In its previously studied form, the method had a trivial zerothorder predictor. In this paper, a more sophisticated firstorder predictor was considered. The cost of such predictor is solution of one linear system per continuation steps.
The firstorder predictor was compared to the zerothorder one on two test cases with finite volume and mimetic finite difference discretizations. Both test cases featured highly nonlinear constitutive relationships. In most cases, the firstorder predictor reduced computational time by decreasing number of continuation steps and Newton iterations.
A part of the work was testing the mimetic finite difference discretization scheme in terms of nonlinear solvers performance. Surprisingly, the scheme was the fastest to produce solution for the capillary barrier test. This fact is attributed to small number of continuation steps. For the realistic test case, it was the only scheme besides TPFA which completed solution. For the capillary barrier test the firstorder predictor gave no advantage, while for the realistic test case it reduced computational time significantly.
Acknowledgements
The reported study was partially funded by Russian Foundation for Basic Research (RFBR), project number 203190126.
References
 (1) L. A. Richards, Capillary conduction of liquids through porous mediums, Physics 1 (5) (1931) 318–333.
 (2) J. Bear, A. H.D. Cheng, Modeling groundwater flow and contaminant transport, Vol. 23, Springer Science & Business Media, 2010.
 (3) M. A. Celia, E. T. Bouloutas, R. L. Zarba, A general massconservative numerical solution for the unsaturated flow equation, Water Resources Research 26 (7) (1990) 1483–1496.
 (4) C. Paniconi, M. Putti, A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems, Water Resources Research 30 (12) (1994) 3357–3374.
 (5) M. W. Farthing, C. E. Kees, T. S. Coffey, C. T. Kelley, C. T. Miller, Efficient steadystate solution techniques for variably saturated groundwater flow, Advances in Water Resources 26 (8) (2003) 833–849.
 (6) C. T. Miller, G. A. Williams, C. T. Kelley, M. D. Tocci, Robust solution of richards’ equation for nonuniform porous media, Water Resources Research 34 (10) (1998) 2599–2610.
 (7) J. E. Jones, C. S. Woodward, Newton–Krylovmultigrid solvers for largescale, highly heterogeneous, variably saturated flow problems, Advances in Water Resources 24 (7) (2001) 763–774.
 (8) H.J. Diersch, P. Perrochet, On the primary variable switching technique for simulating unsaturated–saturated flows, Advances in Water Resources 23 (3) (1999) 271–301.
 (9) D. V. Anuprienko, I. V. Kapyrin, Modeling groundwater flow in unconfined conditions: numerical model and solvers’ efficiency, Lobachevskii Journal of Mathematics 39 (7) (2018) 867–873.
 (10) W. Jäger, J. Kačur, Solution of porous medium type systems by linear approximation schemes, Numerische Mathematik 60 (1) (1991) 407–427.
 (11) I. S. Pop, F. Radu, P. Knabner, Mixed finite elements for the Richards’ equation: linearization procedure, Journal of Computational and Applied Mathematics 168 (12) (2004) 365–373.
 (12) M. W. Farthing, F. L. Ogden, Numerical solution of Richards’ equation: a review of advances and challenges, Soil Science Society of America Journal 81 (6) (2017) 1257–1269.
 (13) D. Anuprienko, I. Kapyrin, Nonlinearity continuation method for steadystate groundwater flow modeling in variably saturated conditions, Journal of Computational and Applied Mathematics (2021) 113502.
 (14) D. Anuprienko, Comparison of nonlinear solvers within continuation method for steadystate variably saturated groundwater flow modeling, Russian Journal of Numerical Analysis and Mathematical Modeling 36 (4) (2021) 183–195.
 (15) E. L. Allgower, K. Georg, Introduction to numerical continuation methods, SIAM, 2003.
 (16) M. T. Van Genuchten, A closedform equation for predicting the hydraulic conductivity of unsaturated soils 1, Soil Science Society of America Journal 44 (5) (1980) 892–898.
 (17) Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resources Research 12 (3) (1976) 513–522.
 (18) A. V. Plenkin, A. Y. Chernyshenko, V. N. Chugunov, I. V. Kapyrin, Adaptive unstructured mesh generation methods for hydrogeological problems, Numerical Methods and Programming 16 (2015) 518–533.
 (19) R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, Handbook of numerical analysis 7 (2000) 713–1018.
 (20) I. Aavatsmark, T. Barkve, O. Bøe, T. Mannseth, Discretization on unstructured grids for inhomogeneous, anisotropic media. Part I: Derivation of the methods, SIAM Journal on Scientific Computing 19 (5) (1998) 1700–1716.
 (21) L. Agelas, D. Di Pietro, I. Kapyrin, R. Masson, The mpfa g scheme for heterogeneous anisotropic diffusion problems on general meshes, in: ECMOR XI11th European Conference on the Mathematics of Oil Recovery, European Association of Geoscientists & Engineers, 2008, pp. cp–62.
 (22) L. Agélas, D. A. Di Pietro, J. Droniou, The g method for heterogeneous anisotropic diffusion on general meshes, ESAIM: Mathematical Modelling and Numerical Analysis 44 (4) (2010) 597–625.
 (23) C. Le Potier, Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés, Comptes Rendus Mathematique 341 (12) (2005) 787–792.
 (24) I. Kapyrin, A family of monotone methods for the numerical solution of threedimensional diffusion problems on unstructured tetrahedral meshes, in: Doklady Mathematics, Vol. 76, Nauka/Interperiodica, 2007, pp. 734–738.
 (25) A. Danilov, Y. V. Vassilevski, A monotone nonlinear finite volume method for diffusion equations on conformal polyhedral meshes.
 (26) O. Misiats, K. Lipnikov, Secondorder accurate monotone finite volume scheme for richards’ equation, Journal of Computational Physics 239 (2013) 123–137.
 (27) K. Nikitin, K. Novikov, Y. Vassilevski, Nonlinear finite volume method with discrete maximum principle for the twophase flow model, Lobachevskii Journal of Mathematics 37 (5) (2016) 570–581.
 (28) K. M. Terekhov, B. T. Mallison, H. A. Tchelepi, Cellcentered nonlinear finitevolume methods for the heterogeneous anisotropic diffusion problem, Journal of Computational Physics 330 (2017) 245–267.
 (29) K. Lipnikov, G. Manzini, M. Shashkov, Mimetic finite difference method, Journal of Computational Physics 257 (2014) 1163–1227.
 (30) L. B. da Veiga, K. Lipnikov, G. Manzini, The mimetic finite difference method for elliptic problems, Vol. 11, Springer, 2014.
 (31) K. Lipnikov, D. Moulton, D. Svyatskiy, New preconditioning strategy for Jacobianfree solvers for variably saturated flows with Richards’ equation, Advances in Water Resources 94 (2016) 11–22.
 (32) E. T. Coon, J. D. Moulton, E. Kikinzon, M. Berndt, G. Manzini, R. Garimella, K. Lipnikov, S. L. Painter, Coupling surface flow and subsurface flow in complex soil structures using mimetic finite differences, Advances in Water Resources 144 (2020) 103701.
 (33) A. S. Abushaikha, K. M. Terekhov, A fully implicit mimetic finite difference scheme for general purpose subsurface reservoir simulation with full tensor permeability, Journal of Computational Physics 406 (2020) 109194.

(34)
K. Lipnikov, G. Manzini, Discretization of mixed formulations of elliptic problems on polyhedral meshes, in: Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, Springer, 2016, pp. 311–342.
 (35) K. Lipnikov, G. Manzini, J. D. Moulton, M. Shashkov, The mimetic finite difference method for elliptic and parabolic problems with a staggered discretization of diffusion coefficient, Journal of Computational Physics 305 (2016) 111–126.
 (36) L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pacific Journal of Mathematics 16 (1) (1966) 1–3.
 (37) D. A. Knoll, D. E. Keyes, Jacobianfree newton–krylov methods: a survey of approaches and applications, Journal of Computational Physics 193 (2) (2004) 357–397.
 (38) I. Kapyrin, V. Ivanov, G. Kopytov, S. Utkin, Integral code GeRa for RAW disposal safety validation, Gornyi Zhurnal (10) (2015) 44–50.
 (39) GeRa software website, http://en.ibrae.ac.ru/code/gera/, accessed: 20210414.
 (40) Y. Vassilevski, K. Terekhov, K. Nikitin, I. Kapyrin, Parallel Finite Volume Computation on General Meshes, Springer Nature, 2020.
 (41) INMOST, a toolkit for distributed mathematical modeling, http://www.inmost.org/, accessed: 20210414.
 (42) I. Konshin, I. Kaporin, K. Nikitin, Y. Vassilevski, Parallel linear systems solution for multiphase flow problems in the INMOST framework, in: Russian Supercomputing Days, 2015, pp. 96–103.
 (43) N. Li, Y. Saad, E. Chow, Crout versions of ILU for general sparse matrices, SIAM Journal on Scientific Computing 25 (2) (2003) 716–728.
 (44) I. S. Duff, J. Koster, The design and use of algorithms for permuting large entries to the diagonal of sparse matrices, SIAM Journal on Matrix Analysis and Applications 20 (4) (1999) 889–901.
 (45) C. M. Oldenburg, K. Pruess, On numerical modeling of capillary barriers, Water Resources Research 29 (4) (1993) 1045–1056.
 (46) S. W. Webb, Generalization of Ross’ tilted capillary barrier diversion formula for different twophase characteristic curves, Water Resources Research 33 (8) (1997) 1855–1859.
 (47) H.J. G. Diersch, FEFLOW: finite element modeling of flow, mass and heat transport in porous and fractured media, Springer Science & Business Media, 2013.
 (48) I. Kapyrin, V. Suskin, A. Rastorguev, K. Nikitin, Verification of unsaturated flow and transport in vadose zone models using the computational code gera, Vopr. At. Nauki Tekh., Ser.: Mat. Model. Fiz. Protsess (1) (2017) 60–75.