ALGORITHMS FOR SOLVING DARCIAN FLOW IN STRUCTURED POROUS MEDIA

This paper presents several algorithms that were implemented in DRUtES [1], a new open source project. DRUtES is a finite element solver for coupled nonlinear parabolic problems, namely the Richards equation with the dual porosity approach (modeling the flow of liquids in a porous medium). Mass balance consistency is crucial in any hydrological balance and contaminant transportation evaluations. An incorrect approximation of the mass term greatly depreciates the results that are obtained. An algorithm for automatic time step selection is presented, as the proper time step length is crucial for achieving accuracy of the Euler time integration method. Various problems arise with poor conditioning of the Richards equation: the computational domain is clustered into subregions separated by a wetting front, and the nonlinear constitutive functions cover a high range of values, while a very simple diagonal preconditioning method greatly improves the matrix properties. The results are presented here, together with an analysis.


Introduction
It is important to be able to predict fluid movement in an unsaturated/saturated zone in many fields, ranging from agriculture, via hydrology, to technical applications of dangerous waste disposal in deep rock formations.
The mathematical model of unsaturated flow was originally published by L. A. Richards [2].The Richards equation problem has undergone various investigations and numerical treatments.Its finite element solution was originally published by Neuman in 1970 for several engineering applications, e.g.dam seepage modeling, see [3].The existence and the uniqueness of its solution was discovered 10 years later, by Alt & Luckhaus [4].A fundamental work analyzing a mass conservation numerical method for the Richards equation was published in 1990 by Celia et al. [5].An analysis of the Richards equation problem has been the subject of several works by Kačur, e.g.[6].Various methods for numerical treatment of the problem have recently been published, see e.g.Kees et al. [7], Šolín et al. [8] or Kuráž et al. [9,10].
An efficient algorithm for numerical solution of the Richards equation should handle temporal adaptivity.This topic has already been discussed by a variety of authors.The first work maintaining the truncation error of the mixed form of the Richards equation was published by Kavetski et al. [14].Tocci et al. [15] proved that the h-based Richards' equation with proper temporal adaptivity can provide a mass conservative numerical solution, and later Miller et al. [16] published an algorithm for temporal adaptation of the step and degree of approximation.
This paper is organized as follows.The first part discusses problems related to inaccurate time integration of the Richards equation.A new and efficient automated time step selection method is presented here, and its efficiency is evaluated on a simple benchmark problem.The second part of the paper discusses poor conditioning of the Richards equation problem.A simple and extremely efficient matrix preconditioner based on diagonal scaling greatly improves the conditioning, and is presented here.Diagonal preconditioning, or more precisely Jacobi preconditioning, is one of the first preconditioning methods that ever appeared in computational history, and its theoretical background was originally published more a hundred years ago.We have not succeeded in our search for the original reference.This method is a very classical approach, and it is even a subject of early courses in applied mathematics.However, it has been found that this standard method is highly efficient for solving the systems of linear equations that arise from numerical approximation of the Richards equation.It has been observed that this preconditioning method even enables very fast convergence of the conjugate gradient method for normal equations, despite powering the condition number.Several more or less recent remarks on using the diagonal preconditioner for CFD and porous media transport problems can be found in [20] and [21],but we have not succeeded in a search for a publication demonstrating the high efficiency of this algorithm as observed.The efficiency of this method will be presented on the basis of one theoretical application and one practical application.However, during the last two years we have not discovered any Richards equation problem where this preconditioning method performed in a significantly different way than as presented here.

Mathematical Model
The problem of Darcian flow in a structured porous medium is usually expressed by the dual porosity conceptual approach.The governing equations for variably saturated Darcian flow and contaminant transport in the dual flow regime were originally published by Gerke, et al. [11].
The following section introduces the strong and weak formulation of dual porosity variably saturated flow.A very brief description is given of the constitutive relations and coefficients involved in this problem, and an interested reader can find more information in the references provided below.

Strong and Weak Formulations
The mathematical problem to be solved is the Richards equation with the dual porosity conceptual approach.
Let Ω be a bounded domain in n (n = 1, 2), with a Lipschitz boundary Let the same hold for the pair Γ D f , Γ N f .Let n denotes the outer unit normal to ∂Ω.For positive The Richards equation in the dual regime was presumed by Gerke & Genuchten [11] as where the subscripts f and m denote the subsystem of fractures (macropores) and matrix blocks (micropores), h is the capillary pressure head function , where θ(h) is the water content function [-], S s is the specific aquifer storage [L −1 ], θ S is the saturated water content [-], ω f denotes the volume ratio of the fracture pore system related to the total pore system, and α w is the first order mass transfer coefficient [L −1 T −1 ] presumed as where β is the dimensionless geometry coefficient, α DP is the characteristic half width [L] of the matrix block, K a is the effective hydraulic conductivity [L T −1 ] of the matrix at or near the fracture/matrix interface, and γ w is the dimensionless scaling factor.The constitutive relation for function θ(h) was supplied by the van Genuchten law [12], and for the K(h) function it was supplied by the Mualem law [13], see Figure 1 for a plot of these functions.
The initial conditions are stated as and the boundary conditions We assume that all coefficients are sufficiently smooth. Let Let the linear space V 1 be a closure of E 1 (Ω), and let V 2 be a closure of E 2 (Ω), in the norm of W 1,2 (Ω). Functions m × [0, T ) and similarly for h f , constitute a variational solution of system (1) together with initial condition and boundary conditions (3)- (4) iff and hold for every The time derivative in equations ( 6)-( 7) was handled with the backward time difference.The continuous unknown h n+1 m,f will be replaced by a piecewisepolynomial function (in our implementation the degree of the polynomial is p = 1) and the test functions ϕ 1 and ϕ 2 with finite element test functions ϕ i 1 and ϕ i 2 to obtain for all i = 1; 2; . . .; D m , D m is the dimension of the finite element subspace of V 1 , and for all i = 1; 2; . . .; D f , D f is the dimension of the finite element subspace of V 2 .The superscript n enumerates the time steps, and ∆t is the length of the time step.
The nonlinear operator was treated by the standard Picard method.This method performs a simple linearization of nonlinear product terms by moving all of them but one to the previous iteration level.
For equations ( 9)-( 10) the following set is obtained (11) for all i = 1; 2; . . .; D m , and for all i = 1; 2; . . .; D f .Here the subscript k denotes the iteration level.The stopping criterion was obtained by looking at the distance of two consecutive Then the only acceptable solution in the next time level is a solution that does not exceed the ranges defined by some critical error of the first order Taylor series approximation of solution h at the previous time level.
Picard approximations of the solutions in the both domains - where ε is the desired accuracy of the method.Equations ( 11)-( 12) constitute the system of linear algebraic equations to be solved in each iteration.The solution of this problem is scattered into h m , and h f .However the scattered outputs h m , h f are not as interesting as the behavior of the capillary pressure head h, which is obtained by averaging the matrix domain m and the fast domain f : One could also hardly find any real problem, where such scattering could be required in the boundary and initial condition setup.Thus the initial condition and the Dirichlet boundary condition are always considered identical for the both domains, and the Neumann boundary condition (volume flux) is always proportionally scattered into the both domains.

Time step adaptivity
It is well known that incorrect discretization of the mass term (time derivative term) yields mass balance errors.This topic, particularly for the Richards equation, was already documented by Celia, et.al. [5] in the 1990s.The primary solution of the Richards equation in h-based form is the capillary pressure head h.The Richards equation law originates from the law of mass balance, and thus the time derivative term can be referred to as a mass term.The original form of the Richards equation published by Richards [2] considers the mass term to be of the form ∂θ(h) ∂t , where θ(h) is the water content function (mass function).

The so-called h-based form of the Richards equation considers the mass term in the form C(h) ∂h ∂t , where C(h) = θ(h) (if zero specific storage is assumed).
Due to the non-linearities involved in the constitutive mass function, see Figure 1, an efficient and accurate time integration algorithm based on the Euler method should handle the selection of an adaptive time step.
Time step adaptivity has been the subject of several works.Most of them focus on limiting the truncation error in the mass term.The original work presenting an implementation of the method of lines for limiting the truncation error of the mixed form of the Richards' equation was published by Kavetski et al. [14].Based on Celia's work [5] the h-based form of the Richards' equation, if numerically approximated by the implicit Euler method, is inappropriate due to the huge inaccuracies in the mass balance.However, it was later presented by Tocci et al. [15], that limiting the truncation error in function h by adaptive time step selection preserves the mass balance and the h-based form of the Richards equation becomes a suitable governing equation for subsurface transport processes.In the works of both Kavetski et al. and Tocci et al., the truncation error was limited by applying the method of lines with a fixed order of approximation.Even more sophisticated work was later published by Miller et al. [16], where temporal adaption was accomplished by variable order and variable step size approximations.
The proposed time integration method does not limit the truncation error of the time derivative term, as mentioned in the previous paragraph.It limits the error of the h-based approximation in a discrete analogue of the mass term ( ∂θ (h)  ∂t ≈ ∆θ(h) ∆t ≈ C(h) ∆h ∆t ).This error can be controlled by defining the ranges for the shift of the solution h in two sequential time steps.The method is depicted in Figure 2. The range of the accepted shift of the solution h in two sequential time steps is a function of h.It will be referred to as h high and h low , and thus the following implicit functions are formed: where ε is the error constant -an error of the first order Taylor series approximation of the constitutive mass function, see Figure 2.This implicit equation was solved by the bisection method.Because it is a resource-consuming computer operation, the equation is tabelized by the initialization procedures.Starting with some initial time step guess, the solution is checked always at the first and last Picard iteration level (the last iteration level of the Picard method is the level, where the distance between two consecutive Picard's approximations is lower than a certain value ( h n+1 k+1 − h n+1 k 2 < ε -see section 2.1).The range h low , h high is evaluated for all mesh nodal values by solving equation (15), where value h is supplied by the previous time step function.It is apparent that evaluating this implicit function directly at each time step for each nodal value would markedly influence the performance of the code, and therefore an efficient and simple algorithm that searches the values in a table was constructed.
The solutions of (15) (h low and h high ) are precalculated for several discrete values of the h function within some user defined ranges and density.Then the solution of ( 15) for any arbitrary h value is obtained by the following algorithm: • the h value is converted into an integer value that points the position of the nearest smallest discrete value of the tabelized function (15), • then the values of solutions h low and h high are obtained by a linear interpolation between the nearest discrete values of the tabelized function (15).
In order to speed up the code performance this technique is also offered to evaluate values for all material functions with real exponents.The user can switch between this approximate and direct material function evaluation method.This method was tested on a standard laptop (Intel CPU 1.6 GHz) for a mesh of 10 6 elements with the nine points integration method.If constitutive functions were evaluated directly, the finite element method matrix was built in approximately 20 min.; if constitutive functions were tabelized as described, the solution matrix was created in approximately one minute.The array initialization took not more than a few seconds.This experiment involved only the van Genuchten's constitutive functions without evaluating the function (15) [24].
If some nodal value in a newly evaluated solution falls beyond the range h low , h high , the new time step solution is rejected, the time step is decreased, and the new solution is tested again -until the criterion is fulfilled.It can be assumed that the constant ε controls the error caused by discretization of the mass term in the h-based form of the Richards equation.The method will be further referred as to the Mass Curve Zone Approach.
It should be emphasized that equation ( 15) has a unique solution only for monotone convex retention curves, e.g. the Gardner's model [17].For the van Genuchten function this algorithm uses only the solution h low , resp.h high , which is closer to the particular evaluated h value.

Case study
A simulation of a vertical one-dimensional infiltration experiment was conducted, see Figure 3.For this numerical example, the dual porosity approach was found to be useless for demonstration purposes.The problem was therefore evaluated on the classical Richards equation only.The selected flow domain Ω was a porous material block 5 m in length and of infinite width.Thus only vertical one-dimensional flow was taking place.Porous material properties (parameters of the van Genuchten and Mualem constitutive functions) were assumed for highly permeable sandstone (α = 0.075 cm −1 , n = 1.89, m = 0.47, K s = 4.42 cm h −1 ).
The domain was discretized by a uniform mesh of grid size 5 cm.
The non-linear Richards equation operator was linearized by the Picard method, the stopping criterion was defined as where subscript k is an iteration level.
In order to evaluate the method presented here, the simulation was performed three times, each time with a different temporal adaptivity mechanism.All selected adaptivity methods were purely empirical. (1.) Number of iterations of the Picard method (max. 10 iterations with a defined stopping criterion)h-based form of the Richards equation.
(2.) Number of iterations of the Picard method (max.10 iterations with the defined stopping criterion)mixed form of the Richards equation (the discrete analogue of the mass term was considered as ∆θ(h) ∆t ). (3.) Mass Curve Zone Approach (the maximal error of a first-order Taylor series, see Figure 2, was 5•10 −6 ).
The initial condition was defined as The top boundary condition was defined as And the bottom boundary was the homogeneous Neumann boundary condition (in technical practice this type of the bottom boundary setup is referred to as free drainage).

Results and Conclusions
As stated in [5] the successful iteration criterion for adaptive time-step selection of the h-based Richards equation is not a sufficient condition to ensure accuracy.The mixed form of the Richards equation provides very accurate results, but the computational effort is excessive in the case evaluated here.The Mass Curve Zone Approach provides results of similar quality, but the total number of iterations (number of calls of the linear equation solver) is only approximately 4 % of the amount required for solving the mixed-form Richards equation.This method offers an efficient condition for a numerical solution of the h-based equation based on the implicit backward Euler scheme, exactly reflecting the problem of the h-based approximation.Compared to the Tocci approach [15], our algorithm does not limit the truncation error of the time derivative term, it only limits the error of the discrete analogue of the mass term if the h-based form of the Richards' equation is considered.However, our algorithm behaves in a different way.It should be evaluated whether there is any relation between these two errors, or at least which criterion is more appropriate for achieving mass consistency.
For particular results see Table 1.

Preconditioner
Finite element approximation of the Richards equation problem combined with the fully implicit Roth method typically forms very badly-conditioned matrices, especially if coarse-grain material is involved.The plot of the unsaturated hydraulic conductivity, Figure 1 right, covers a typical capillary pressure head range for a certain Richards equation problem, and thus it is not uncommon to solve problems where the non-linear model parameters can vary between eight or even ten orders of magnitude.The poor conditioning of the Richards equation is caused by the huge range of values of the constitutive functions.Another property of the Richards equation also applies.Due to the small time steps that are typically required for an accurate solution of this problem, and due to lumping of the capacity matrices in order to preserve the discrete maximum principle, the matrices of its numerical approximation are diagonally dominant.Then it was observed that badly conditioned matrices of this type are formed out of clusters with a lower condition number.For simple infiltration problems, the clusters are formed at (1.) fracture and matrix flow domains, (2.) regions separated by the wetting front.
And thus it can be summarized that, • due to the small time steps, matrices are diagonally dominant, • due to the small time steps and the continuity of the solution (it is in fact just a non-linear diffusion with convection), the matrices are "similar" in two sequential time steps, • the domain can be clustered into subdomains with similar values of the constitutive functions, and thus the matrix can be clustered into "submatrices" with a lower condition number (as the problem is nonlinear and often unsteady, the position of these subdomains is not constant in time).
Based on these assumptions, the perfect preconditioner should adjust the material properties in the computational domain into such a state, that the material parameters (values of the constitutive functions K(h, x) and C(h, x))are transformed into certain homogeneous values.Behavior of this kind is actually achieved by the standard Jacobi preconditioner (diagonal scaling) assumed from the right hand side.The Jacobi preconditioner is one of the simplest forms of preconditioning, in which the preconditioner matrix is the diagonal of the original matrix A. The method can be described by the following scheme . . .1.0 and thus this method can be referred to as a right hand side diagonal preconditioner where This method suffers from one notable disadvantage -it creates non-symmetric matrices even out of symmetric matrices.However, transport processes are in general a convection-diffusion process, and thus they are non-symmetric.In the context of the conjugate gradient method for normal equations, this method of preconditioning can be formulated as D T A T AD, which is actually a standard approach.

On Selecting an Algorithm for Solving the Systems of Linear Equations
Transport processes in porous media are in general of a convection-diffusion-reaction form.In the case of the standard Richards equation the convection is caused by gravitational forces, and then the problem can be formulated symmetrically with respect to the total hydraulic head H instead of the pressure head h.This approach would be efficient for the case studies evaluated here, but the transport processes in porous media are more accurately described by a coupled process of moisture, heat and solute transport, as described by Noborio et al. [18].Thus the search for an efficient solver should not be limited to symmetric problems only.Moreover, preserving the formulation of the Richards equation in the pressure head form h enables density dependent transport problems to be solved, e.g.sea water intrusion.This form was therefore preferred in implementing the DRUtES code.Secondly, the diagonal scaling method that is applied forms non-symmetric matrices out of symmetric matrices.This property is obviously inconvenient, but based on the following assumption the cost for losing the symmetricity, either due to a less appropriate formulation of the Richards equation or due to the application of right hand side diagonal scaling, is insignificant.
Let A be the matrix obtained by a numerical approximation of the Richards equation problem into the system of linear equations Ax = b, as described in section 2.1.For each iteration of the Picard method, matrix A and vector b is assembled from equations ( 11)-( 12) and from the boundary condition setup (4).Due to nonlinearities in the Richards equation, and due to preserving the mass balance, the problem is temporally discretized into small time steps.Together with lumping the capacity matrix as recommended by Rank et al. [19], in order to preserve the discrete maximum principle, we could start with the hypothesis of a strong diagonal dominance of matrix A.
Then matrix A can be decomposed into the sum of its diagonal and nondiagonal part where D = diag(A), and obviously B = A − D. The right hand side diagonal scaling rearranges the original problem Ax = b into where x = D −1 y.Then the matrix AD −1 can be considered as where I is an identity matrix.If we apply the conjugate gradient algorithm for normal equations to the matrix I + BD −1 from ( 24), then the following is obtained where G is the matrix that is processed in the conjugate gradient algorithm.
Based on the hypothesis of the strong diagonal dominance of matrix A, the Frobenius norm of and thus it could be concluded that Frobenius norm of the matrix G from (25) must be close to 1.0.This conclusion will be deeply explained.Let us consider matrix F as and thus Then the Frobenius norm of matrix F is referred to We have already assumed the Frobenius norm of BD −1 to be small, see (26), due to the triangle inequality and due the submultiplicativity of the Frobenius norm we can conclude that In order to estimate the efficiency of the conjugate gradient solver applied to the matrix G from (25), we should estimate the condition numbers of matrix F from (28).It is well known that the absolute values of all eigenvalues are lower than its matrix norm (this applies for all commonly used norms including the Frobenius norm).And thus the following applies And thus matrix G from (25), which is related to matrix F by definition (27), has eigenvalues within the range The condition number is evaluated as where κ is the condition number of matrix G, λ F min is the minimal eigen value of matrix F, and λ F max is the maximal eigen value of matrix F, and thus In order to fulfill the proposed prerequisite, the norm BD −1 must be small, but as already mentioned the value of this norm for our Richards equation problems was typically observed within the range The conjugate gradient method for normal equations is not a typical choice due to squaring the condition number.The typical solution is based on GMRES or BiCGSTAB.GMRES is an optimal method on the Krylov subspace generated by the initial residual r 0 and matrix A. But its great disadvantage is that long recurrences are required, so it is almost impossible to use the original form of GMRES, as the long recurrences may cause serious problems with computer memory allocation.GMRES is typically used in a variety of modifications, as restarts or shorten recurrences.As presented by Greenbaum et al. [22] the convergence rate does not depend on the eigenvalues, and thus achieving an accurate solution in restarted GMRES is uncertain (at least on theoretical problems).The BiCGSTAB method uses short recurrences, but its great disadvantage is the possibility of a so-called break-down.
The simple conjugate gradient method for normal equations does not suffer from any these disadvantages.The only disadvantage is squaring the condition number, but on the basis of the given hypothesis for this particular class of problems worsening the condition number is negligible, and thus this only disadvantage is absent here.
The conjugate gradient method for normal equations is implementary simple, and as explained for diagonally scaled diagonally dominant problems, it is also an efficient algorithm.This method was implemented into DRUtES code, and bas been used successfully for a variety of theoretical and practical problems.

Theoretical Case Study
This example is a purely theoretical case study.The matrix and fracture domain parameters are consid- ered identical here, and thus the problem is described by the classical Richards equation model only.The material parameters are depicted in Table 3, 3rd line "gravel bar.".The boundary condition setup is stated as: The initial condition setup is stated as: The domain scheme is depicted in Figure 4.
The spatial discretization was a uniform triangular mesh with 10201 nodes and 20000 elements.

Conditioning Improvements
The diagonal scaling algorithm greatly improved the problem conditioning.The overall number of Picard iterations -linear equation solver calls -was 49001.The matrix of dimension 10100 was solved by the conjugate gradient method for a normal equation with an average of 33 iterations, and the maximal number of iterations required per CG solver call was 67.The stopping criterion was the residual size r k < 10 −12 .Table 2 depicts the results that were obtained.Figure 5 depicts the plot of the conjugate gradient iterations for preconditioned and non-preconditioned matrices.As expected, the CG method for a normal equation without preconditioning converges so slowly that its application can be concluded as an inappropriateconsidered inappropriate.This is a wellknown fact, but simple diagonal scaling enables fast convergence.in the Bídnice hill, 40 m below the surface and 70 m above the ground water table, and it is surrounded by tectonically fractured claystone.Under normal circumstances, this is a very stable unsaturated environment.This is therefore a typical problem for applying the Richards equation model.In order to protect its surroundings from radiotoxic pollution, the repository is shielded by a 30 cm thick gravel layer that acts as a capillary barrier.The capillary barrier effect is based on the unsaturated hydraulic properties of gravel.The unsaturated hydraulic conductivity of gravel reaches very high values at and near full saturation, but with decreasing saturation it falls steeply to almost zero values.

Practical Case study
The problem was conceptualized by model (1), where the tectonic fractures are represented by the fast domain (f ), and the rock matrix is represented by the matrix domain (m), the material properties are depicted in Table 3, and the domain scheme is depicted in Figure 6.The coupling parameters of the model (1) were identified in [9], and are stated as: α = 8.89911 cm −1 d −1 , and ω f = 0.0292884.The hydraulic parameters (parameters of the van Genuchten and the Mualem function) are depicted in Table 3, and Figure 1 is a plot of the constitutive functions.The repository interior -concrete, and the repository barrier -gravel, were obviously formed out of identical material in both the matrix domain and the fracture domain.The different material conditions scattered into the matrix and fracture domains appeared in the rock material.
The initial condition was defined as The domain was discretized by a non-uniform triangular mesh that consists of 12027 nodes and 23771 elements.

Numerical Properties
The material properties of this problem cover several orders of magnitude, and thus finite element approximation forms ill-posed matrices.In order to evaluate the efficiency of the preconditioner that is applied here, several early states were solved both with and without preconditioned matrices.When the stopping criterion definition from the previous paragraph was used, it was observed that for non-preconditioned matrices the number of iterations reached an average   value of 17687.In the case of a preconditioned state, the average number of iterations was 18.The matrix dimension was 23898.Thus, for this particular case, this simple preconditioner decreased the total number of iterations required to solve the unsteady problem by a multiple of approximately 10 −3 .As was mentioned above, this proposed method suffers from a noticeable disadvantage -for symmetric problems the matrix is converted into a non-symmetric matrix.Thus the diagonal scaling was also assumed in a different way -from the left hand side and then the scaling can be described by the following scheme It is well known, see eg. [23], that if A were a symmetric positive definite matrix, then the left and right hand side preconditioners would be equivalent.It should be further noted, that the left hand side preconditioner in the context of the conjugate gradient method for normal equations forms A T D T DA, and this form is not commonly used.
This preconditioner was again evaluated in our repository study.It was observed that the left hand and right hand side preconditioners both behave similarly if the solution is "smooth", while for less smooth solutions the right hand side preconditioner offers better performance, see Figure 7.The "smoothness" of the solution was defined by the criterion where n is matrix dimension, x is vector of solution (Ax = b), and x i is some particular value at position i from the solution vector x.Constant s m is only a certain indicator of the "smoothness" of the values in vector x.Then, as observed (see Figure 7), at time e.g.t = 2 days → s m ≈ 10 4 -the behavior of the two methods was comparable.However, for some later simulation time levels the right hand side preconditioner offered better performance, e.g. at time t = 80 days → s m = 10 6 .If the right hand side preconditioner is used, together with the small time steps involved, the "smoothness" of the solution remains hidden in the diagonal matrix D. This is because the solution x prcd of matrix A prcd is related to solution x of matrix A by and thus for the small time steps involved, together with an assumption that the solution of the Richards equation is continuous, the s m constant from equation (42) of the x prcd vector reaches a low value.

Conditioning Improvements
The unsteady problem of infiltration around the repository can be divided into the following periods (1.) intense initial infiltration -the Dirichlet boundary condition does not satisfy the initial condition.Such discontinuity causes intense processes in the initial states, see Figure 8; (2.) the wetting front is getting closer to the capillary barrier shield, the repository is being bypassed, see Figure 9; (3.) the wetting front has reached the repository shield, steep gradients appear at the wetting front/barrier interface, see Figure 10.
The right hand side preconditioner efficiency was evaluated for each period, for a selected time and iteration level.Minimal and maximal eigenvalues were evaluated for the matrix A T × A, as the problem is non-symmetric, and thus a conjugate gradient algorithm for normal equations was applied.Maximal eigenvalues were evaluated by the power method and minimal eigenvalues were evaluated by inverse iteration with Cholesky decomposition.
The results are depicted in Table 4.

Conclusions
This paper has aimed to present several new algorithms that were implemented in the DRUtES opensource project [1] -a finite element solver for nonlinear coupled parabolic problems, namely the Richards equation with the dual porosity approach, in order to improve the accuracy and efficiency of its numerical treatment.
The Mass Curve Zone Approach, an efficient and accurate method for the Euler time integration method, has been presented and tested.This method is based on an automatic time step selection.The length of the time step is selected in such a way that the error of the discrete approximation of the time derivative term (mass term) is always kept below a certain value.
A different problem concerned the conditioning of the Richards equation problem.It is well known that finite element approximation of the Richards equation forms badly-conditioned matrices.A very simple and efficient preconditioning method based on diagonal scaling has been presented.Then the method of preconditioning was tested on a real problem -infiltration into a capillary barrier layer that was used as a protective shield around a waste repository.The method presented here efficiently decreased the condition number of the evaluated problem by more than 13 orders of magnitude.The best performance was observed for the initial time levels, where the condition number of the original matrix was κ = 1.43 • 10 13 , while after preconditioning κ = 1.00215, and thus the preconditioned matrix was even close to an identity problem.Generally speaking, all matrices that appeared in this simulation experiment had a condition number around 10 13 -10 15 , and the proposed preconditioner decreased the conditioning to values lower than 10 2 -10 3 .The efficiency of the proposed preconditioner is highly dependent on the selected time step.
All simulated cases demonstrate great effectiveness of the discussed automated time step selection method for Euler time integration, and for the proposed matrix preconditioner.Generally speaking, both the design and the implementation of the preconditioner are trivial, and it is therefore a recommended method for conjugate gradient based finite element solvers of transport processes in a porous medium.

Time
No

Figure 1 .
Figure 1.Plot of the constitutive functions for the Richards equation (1).Left -the van Genuchten functionretention curve -relation between water content θ and the capillary pressure head h (mass function).Right -the Mualem function -unsaturated hydraulic conductivity)

Figure 2 .
Figure 2. The critical point of solution h was assumed from nodal values of the discretization mesh.Then the only acceptable solution in the next time level is a solution that does not exceed the ranges defined by some critical error of the first order Taylor series approximation of solution h at the previous time level.

Figure 3 .
Figure 3. Scheme of the one-dimensional vertical infiltration experiment.

Figure 4 .
Figure 4. Scheme of the two-dimensional theoretical experiment.

Figure 6 .
Figure 6.Domain scheme of the repository problem.

Figure 7 .
Figure 7. Number of iterations required by the conjugate gradient algorithm to converge for the repository problem.

Figure 8 .
Figure 8. Intense initial infiltration -the Dirichlet boundary condition does not satisfy the initial condition, such discontinuity causes intense processes in the initial states.

Figure 9 .
Figure 9. Wetting front is getting closer to the capillary barrier shield, the repository is being bypassed.

Figure 10 .
Figure 10.Wetting front has reached the repository shield, steep gradients appears at the wetting front/barrier interface.

Table 1 .
Summary of code execution after the numerical infiltration experiment.

Table 2 .
This example is motivated by an existing nuclear waste repository for low-level and intermediate-level waste.The name of the repository is Richard and it is situated near Litoměřice, Czech Republic.The repository is located within a deep rock formation Conditioning improvements before and after right hand side diagonal scaling.If λmax < 1.0, then most likely more iterations are required to improve the accuracy of the eigenvalue estimation.

Table 3 .
The material properties table.