First-order continuation method for steady-state variably saturated groundwater flow modeling

by   Denis Anuprienko, et al.

Recently, the nonlinearity continuation method has been used to numerically solve boundary value problems for steady-state Richards equation. The method can be considered as a predictor-corrector procedure with the simplest form which has been applied to date having a trivial, zeroth-order predictor. In this article, effect of a more sophisticated predictor technique is examined. Numerical experiments are performed with finite volume and mimetic finite difference discretizations on various problems, including real-life examples.


page 15

page 16


Comparison of nonlinear solvers within continuation method for steady-state variably saturated groundwater flow modeling

Nonlinearity continuation method, applied to boundary value problems for...

Nonlinearity continuation method for steady-state groundwater flow modeling in variably saturated conditions

Application of nonlinearity continuation method to numerical solution of...

A Fast Method for Steady-State Memristor Crossbar Array Circuit Simulation

In this work we propose an effective preconditioning technique to accele...

Mathematical modelling of an unstable bent flow using the selective frequency damping method

The selective frequency damping method was applied to a bent flow. The m...

An explicit predictor/multicorrector time marching with automatic adaptivity for finite-strain elastodynamics

We propose a time-adaptive predictor/multi-corrector method to solve hyp...

Field model for complex ionic fluids: analytical properties and numerical investigation

In this paper, we consider the field model for complex ionic fluids with...

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 steady-state cases are subject of extensive research farthing2003efficient ; farthing2017numerical .

This paper focuses on improvement of a special solution approach for boundary value problems for the steady-state 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 predictor-corrector continuation methods with a trivial, zeroth-order predictor step and the corrector step being application of an iterative nonlinear solver.

This paper contributes to the nonlinearity continuation method for the steady-state Richards equation in the following ways.

  1. The nonlinearity continuation method is formulated in predictor-corrector terms, the first-order predictor is applied and its effect on performance is examined;

  2. 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 predictor-corrector 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


and the Darcy law for water flux


the steady-state Richards equation richards1931capillary ; bear2010modeling


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 steady-state 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 steady-state Richards equation in the following form:


The same applies to the flux expression (2), and it becomes


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 steady-state 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 :


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 non-physical 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 diffusion-type 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 two-point flux approximation (TPFA) remains the most popular choice for its simplicity. However, TPFA does not give approximation for non--orthogonal grids, which are very common in practical applications. Other FV schemes include linear multi-point flux approximations aavatsmark1998discretization ; agelas2008mpfa ; agelas2010g and schemes with nonlinear two- and multi-point 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:

The relative permeability for a face can be either calculated as a half-sum 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


where is a matrix with solution-dependent 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 first-order PDEs (1) and (5). These equations represent diffusion problem in mixed formwith nonlinear solution-dependent coefficient. The scheme used in this work chooses the divergence as the primary operator and uses finite-volume-like 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 solution-dependent 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 cell-centered hydraulic head unknowns (one value per cell) and face-centered flux unknowns (one value per face, representing the normal flux). The respective nonlinear system, therefore, has significantly larger size compared to cell-centered 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 non-convex 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 predictor-corrector 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



is the respective vector of discrete unknowns (for finite volume schemes it consists of cell-based head unknowns, while for the mimetic scheme it also contains face-based 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 steady-state 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


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 Predictor-corrector formulation

The nonlinearity continuation procedure described above lies in the class of predictor-corrector 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 zeroth-order predictor: it simply sets solution of the previous problem as initial guess for the next problem.

Find : with an iterative solver;
while  do
       while  do
             1. Predictor step:
             set initial guess for ;
             2. Corrector step:
             Find : with an iterative solver;
             if solved then
             end if
       end while
      if didn’t find  then
       end if
end while
Algorithm 1 Nonlinearity continuation method in predictor-corrector terms

4.3 First-order predictor

In algorithm 1 the predictor step is setting initial guess for . Previously, a trivial, zeroth-order predictor was used, which has the form


A more sophisticated predictor may be used. In this so-called first-order 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


where is the Jacobian matrix calculated for the last obtained solution, and the right-hand side is approximated with a finite-difference 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 first-order 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 first-order predictor on the number of continuation steps and nonlinear iterations, compared to the zeroth-order 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 ; gera-site . GeRa was originally designed for safety assessment of radioactive waste repositories and nowadays is a general-purpose hydrogeological modelling tool. Computational core of GeRa is based on INMOST vassilevski2020parallel ; inmost-site , a platform which provides tools for numerical modeling on general unstructured grids with built-in MPI parallelization.

Of the discretization schemes used, only TPFA and MPFA-O are available in GeRa. NTPFA-B and NMPFA-B 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 Bi-CGSTAB preconditioned by second order Crout-ILU decomposition konshin2015parallel

with inverse-based 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 two-layered 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 two-dimensional 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.

Figure 1: Computational domain and boundary conditions for the capillary barrier test (stretched 5 times vertically)
, m/day , 1/m
Sand 18.144 3.9 5.74 0.154 0.39
Gravel 8640 490 2.19 0.011 0.42
Table 1: Media paramaters for the capillary barrier test

Numerical experiments were conducted on a 3200-cell 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 first-order predictor was able to decrease total computation time by up to 40% by reducing number of continuation steps and Newton iterations. For the NMPFA-B scheme, the first-order 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 first-order 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 per-iteration 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
MPFA-O, 0 375.0 31 (33) 667
1 281.5 24 (25) 432
NTPFA-B, 0 48.7 35 (37) 409
1 29.3 26 (27) 265
NMPFA-B, 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
Table 2: Comparison of zero- and first-order predictors for the capillary barrier test, power continuation function
Figure 2: Hydraulic head and water saturation distributions for the capillary barrier test (stretched 5 times vertically)

5.3 Realistic case

The second problem involves groundwater flow in a vertical cross-section of a real-life 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 non-smooth 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 non--orthogonality 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 first-order 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 first-order 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.

Figure 3: Computational domain and boundary conditions for the realistic case
, 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
Table 3: Media properties for the realistic case
, 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
Table 4: Comparison of zero- and first-order predictor for the realistic test
Figure 4: Water saturation distributions computed for the realistic case with TPFA (top) and MFD (bottom)
Figure 5: ”Checkerboarding” in water saturation observed for TPFA in a part of the domain

6 Conclusion

The nonlinearity continuation method, which has been recently applied for boundary value problems for the steady-state Richards equation, can be considered as a predictor-corrector procedure. In its previously studied form, the method had a trivial zeroth-order predictor. In this paper, a more sophisticated first-order predictor was considered. The cost of such predictor is solution of one linear system per continuation steps.

The first-order predictor was compared to the zeroth-order 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 first-order 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 first-order predictor gave no advantage, while for the realistic test case it reduced computational time significantly.


The reported study was partially funded by Russian Foundation for Basic Research (RFBR), project number 20-31-90126.


  • (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 mass-conservative 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 steady-state 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–Krylov-multigrid solvers for large-scale, 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 (1-2) (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 steady-state 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 steady-state 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 closed-form 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 XI-11th 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 three-dimensional 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, Second-order 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 two-phase flow model, Lobachevskii Journal of Mathematics 37 (5) (2016) 570–581.
  • (28) K. M. Terekhov, B. T. Mallison, H. A. Tchelepi, Cell-centered nonlinear finite-volume 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 Jacobian-free 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, Jacobian-free 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,, accessed: 2021-04-14.
  • (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,, accessed: 2021-04-14.
  • (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 two-phase 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.