Simulation of flows in heterogeneous porous media of variable saturation

A R T I C L E I N F O A B S T R A C T Article History: Received: 12 June, 2017 Accepted: 15 July, 2017 Online: 26 December, 2017 We develop a resolution of the Richards equation for the porous media of variable saturation by a finite element method. A formulation of interstitial pressure head and volumetric water content is used. A good conservation of the global and local mass is obtained. Some applications in the case of heterogeneous media are presented. These are examples which make it possible to demonstrate the capillary barrier effect.


Introduction
The study of the flow of fluids in porous media of variable saturations is of practical interest in many fields, in agriculture, civil engineering as well as environment. Control of the movement of fluids in industrial waste, tailings and municipal waste sites requires a thorough knowledge of the dynamics of these fluids; either to assess the extent of the damage, or to design retention or containment structures or to implement a rehabilitation strategy.
It is not always technically easy to represent the physical reality of these problems in the laboratory or to design models at a reasonable cost. The mathematical and numerical modeling of the flow of fluids in porous media has for some decades made considerable progress in the wake of the rapid development of the means of calculation. The flow is modeled by nonlinear partial differential equations. These equations rarely have analytical solutions, or at the cost of simplifications which make the model without practical interest. It is rather by the various numerical approximation methods that these equations are solved. The nonlinear character of the properties of the porous media of variable saturation, the high heterogeneity of these media and the mixed nature of certain boundary conditions such as those prevailing on the soil-atmosphere interface, make the problem singularly difficult and require the implementation of robust and efficient approximation methods.
In this study, we propose a mixed formulation in volumetric water content and in interstitial pressure head. The time discretization is performed with an implicit Euler scheme with variable time step. In space, a finite element method is used. The nonlinear problem is solved by the method of successive iterations. Although the method is of the first order, taking into account the variations of the volumetric water content in the iterative scheme makes it possible to obtain the convergence with a good conservation of the mass.
As applications, we first propose the study of the classic problem of infiltrations of water in a very dry environment. It is a question of testing the convergence and the robustness of the proposed method in comparison with the results of the literature. Subsequently, the case of flow in a heterogeneous medium is examined. We study more precisely the flow of water in a medium formed of several layers of materials of different textures. The flow conditions through these layers are often complex due mainly to the non-linear behavior of the various hydraulic variables of the medium and their discontinuity in the passage from one layer to the other. Due to the contrast of the saturation state between the layers, a capillary barrier effect occurs at their interface. This can greatly reduce the movement of fluids (water and air) between layers. This mechanical phenomenon has been extensively exploited in the design of recovery systems for storage sites of various types of waste (industrial, mining, domestic, etc.). These systems contribute to the isolation of these releases from their immediate environment by restricting the movement of water and oxygen. The general characteristics of multi-layer barriers have been dealt with extensively in the mining and geotechnical literature see [1][2][3][4].

Equations of flow
Unsaturated subsurface media are triphase: the medium skeleton is the solid phase, water and air form the liquid and gaseous phases, respectively. When the skeleton is assumed to be rigid and immobile, and air is at atmospheric pressure, the flow of water is governed by Richards' equation [5] . It is a simple (seemingly) equation model describing the evolution of volumetric water content in unsaturated soils. The Richards equation is obtained from the law of conservation of the mass of water and a law of behavior, the law of Darcy [6]. The Richards equation in a mixed formulation, volumetric water content and interstitial pressure head, is: Where Ω is a domain of the plane ℝ 2 (or the space ℝ 3 ) representing (geometrically) the porous medium, θ is the volumetric water content ([L 3 L −3 ] dimensionless), u is the Darcy velocity [L T −1 ], h is the interstitial pressure head [L], f is a source function [T −1 ] and z is the elevation [L] . K is the hydraulic conductivity tensor [L T −1 ] of the medium that is written in the form: Where K A is the anisotropy tensor, K is the hydraulic conductivity, k s is the hydraulic conductivity of the medium at water saturation and k r is the relative conductivity.
To the equation (2.1) are associated constitutive relationships between the volumetric water content and the interstitial pressure head on the one hand, and between the hydraulic conductivity and the interstitial pressure head on the other hand. These are empirical relationships such as those proposed by [7], Van  Where: is the residual water content and is the water content at saturation.

Note 1:
In the literature, the Richards equation is also written and studied in a so-called interstitial pressure formulation. To obtain this formulation, we put: where: C(h) is the hydraulic capacity function of the medium. In order to solve mathematically the partial differential equation (2.1), it is associated with initial and boundary conditions.

Initial conditions
We suppose that: where: h 0 is a given function on Ω.

Boundary conditions
We can associate to the equation (2.1) several types of boundary conditions according to the physical model studied: where: Γ D , Γ N and Γ G are respectively parts of the boundary of the domain Ω and h D , σ n and σ G are given functions and n is the normal external to Γ N and Γ G . In flow models across soils of varying saturation, boundary conditions are not reduced to the previous types. Non-standard conditions are imposed as seepage boundary conditions, atmospheric boundary conditions, and free drainage boundary conditions.
The seepage boundary conditions consist in imposing the interstitial pressure head at atmospheric pressure on parts Γ S , a priori unknown, of a surface, subject to a potential seepage.
The atmospheric boundary conditions are applied to the surface of the soil to take into account the physical phenomena of precipitation and evaporation on the surface Γ . These conditions must respect the constraints of the state of saturation or not of the soil in the vicinity of the surface in the form of the two inequalities: where: E S is the maximum potential flux on the surface and h A is the minimum pressure head under current soil conditions. The boundary conditions of free drainage consist in imposing the Neumann condition with: σ G (x, z, t) = 0.
The finite element method has the advantage of taking into account general boundary conditions and allows to study the equation of the flow in complex geometries. In addition, the finite element method is flexible and practical.
The main difficulty of the numerical resolution of the Richards equation is its nonlinearity in h by the volumetric water content θ(h) and the hydraulic conductivity K(h) on the one hand, and the discontinuities due to heterogeneity of the environment on the other hand. Numerical methods present convergence difficulties as soon as these functions undergo sudden transitions in the vicinity of certain values of the interstitial pressure head. This is the case, for example, during the propagation of the wetting front in dry soil. The solutions exhibit instabilities with bad mass conservation see [9,16,17].
In the literature, the Richards equation is written and discretized under three types of formulations: in h, in θ and in h and θ. These formulations are formally identical, but sometimes produce different numerical results [18][19][20]. The formulation in h does not ensure a good conservation of the mass unless we consider a fine mesh and a small time step. The formulation in θ gives a very good conservation of the mass [9], but it does not make it possible to treat correctly the flow in the saturated zones, because the equation degenerates. The so-called mixed form in h and θ, for its part, ensures a good approximation of the solution and a better conservation of the mass. A discussion of this subject in the one-dimensional case is found in [9], [18]. The choice of iterative methods, such as Newton-Raphson or the method of successive approximations (Picard) as well as their various variants, also depends on the considerations of the cost in computation time.
In this section, we study the approximation of the solutions of the Richards equation. We use a finite difference scheme for time discretization with a variable time step. On each time step, we obtain a nonlinear stationary problem, which is solved after linearization by an iterative method. At each iteration, a linear problem is solved by a standard finite element method.

Time discretization
In the following, we present the time discretization of the Richards equation in the so-called mixed form h and θ: We denote by ∆t the time step and by t 0 , t 1 , … , t n , t n+1 , … the points of the discretization, with t n+1 = t n + ∆t. We also denote by h n , θ n and f n respectively the values of h, θ and f at time t n , ie h n = h(•, t n ) , θ n = θ(h(•, t n )) and f n = f(•, t n ) . The discretization of (3.1) by the implicit Euler scheme leads us to solve the nonlinear problem, finding h n+1 satisfying: To be able to write an iterative scheme of successive approximations of this problem, we must first linearize the function θ(h). If we denote by h n+1,ϑ , h n+1,ϑ+1 the ϑ + 1 th and ϑ th approximation of h n+1 , we then write using the Taylor formula : The linearization of (3.2) is written in the form The iterative process (3.5) is called by some authors a method of modified Picard iterations [9,18].

Note 2:
In the equation (3.2) we could have linearized the function K(h n+1 ), as we did for θ(h n+1 ), by writing its Taylor expansion to the order 1 in the neighborhood of a given approximation h n+1,ϑ , which will naturally involve the derived from dK/dh. This would lead to an iterative second order scheme of the Newton type. From this point of view, scheme (3.5) appears in fact as an approximation of Newton's method. The use of the Newton method is advantageous in the case of very dry soils where it converges faster, but with a higher cost [18].

Finite element discretization
We decompose the domain of the flow Ω into a finite number of elements. The approximate resolution of (3.5) by the finite element method in a standard Galerkin type formulation consists in determining the solution h n+1,ϑ+1 in the form: Where: nn denotes the number of interpolation nodes, {φ 1 , … , φ nn } are the basic functions of the Lagrange interpolation space, and h j n+1,ϑ+1 is the degree of freedom of h i n+1,ϑ+1 at the jth node. Galerkin's method amounts to writing that h i n+1,ϑ+1 verifies (3.5) in the weak sense, that is, the integral identity: for any basic function φ i .
By introducing the explicit form of the terms of identity (3.6), we obtain the linear systems.
We have thus reduced the approximation of the problem (2.1) to the resolution of a linear system of type (3.7) on each time step. The test on θ (step (7)) allows a better control of the volumetric water content and ensures a better convergence and conservation of the mass.

Linear system resolution
The linear system obtained in (3.7) is of the form: Ax = b, where A is a symmetric and definite positive sparse matrix N × N, and b is a vector of ℝ N given, x is the unknown vector of ℝ N . When Matrix A is small, the system can be solved by a direct method such as the Gauss or Choleski elimination method. If, on the other hand, A is large, it is more advantageous to use iterative methods such as the conjugate gradient method or the preconditioned conjugate gradient method. In our case, we use the conjugate gradient method preconditioned by incomplete LU factorization [21].

Numerical simulations
In this section, some examples of flow tests in porous media are presented. The main aim is to demonstrate the effectiveness of the formulation used and the robustness of the finite element method. This is particularly illustrated in the case of water infiltration in very dry porous media as well as in the case of drainage in heterogeneous media.

Infiltration in a homogeneous column
The classical problem of infiltration of water into a homogeneous porous column is considered. This example has been widely used by several authors to validate their numerical methods for solving flow problems in porous media of variable saturation. See [9,18,20,22,23,24]. The hydraulic properties of the medium are:  In this first test, the height of the column is L = 30 cm. The initial condition is h(. ,0) = − 1 000 cm , while the boundary conditions are h(L, t) = − 75 cm at the surface of the column and h(. ,0) = − 1 000 cm on its base. The simulation time is 6 hours. A regular mesh has been produced with the height of the elements ∆z = 0,25 cm. It has been used with several time steps, from the finest ∆t = 0,1 s to the widest 50 s ≤ ∆t ≤ 500 s. The solutions obtained are quite comparable to those of Lehmann in [18], and to the semi-analytical solution of Philip [25] as well as those obtained by [9]. We also considered a fine mesh with elements of size ∆z = 0,1 cm and a time step ∆t = 0,1 s. This mesh was used to calculate the comparison solution. In the case of the fine mesh and the small time step, the same results are obtained with a good conservation of the mass. In figures 3 and 4, the curves of the interstitial pressure head and the curves of the volumetric water content as a function of the elevation for a regular mesh, whereas in figures 5 and 6 the same representations in the case of fine mesh.

Infiltration in a heterogeneous (multilayer) column
An example given by Lehman in [18]. It is a heterogeneous column of height L = 180 cm. In this example, it is sought to test the hydraulic behavior of a heterogeneous medium with boundary conditions relating to a wetting and drainage cycle. The porous medium consists of three layers each 60 cm thick. The hydraulic parameters of the first and the third layer (Berino loamy fine sand) and those of the second layer (Glendale clay loam) are given by: The characteristics of the medium, θ(h) and K(h) are shown in figures 8 and 9.

I) Initially wet medium
The initial interstitial pressure head is taken as h(. ,0) = −100 cm. We first considered a regular mesh, with the height of elements ∆z = 5,0 cm. For the time steps, 0,01 ≤ ∆t ≤ 1,0 was used. Thereafter, a mesh with element height ∆z = 1,0 cm was used. In both cases, we find the same results as Lehman and al. [18], figures 11 and 12. In figure 12, it is noted that the middle layer has been saturated rapidly, while the top and bottom layers desaturate. Thus we are in the presence of an effect of capillary barrier between layers.

II) Initially dry medium
For this test, it is assumed that the three media are initially very dry. The initial interstitial pressure head is taken as h(. ,0) = −10 000 cm and as a time step, successively ∆t = 0,1 days and ∆t = 1,0 days. To obtain good convergence, we have considered the fine mesh (∆z = 1,0 cm) associated with a rather small time step (∆t = 0,0001 days). Figures 13 and 14 shows the interstitial pressure head and the volumetric water content. As in the previous case, a capillary barrier effect is highlighted. Although the medium was initially dry h = −10 000 cm, the layer of fine material in the middle of the column saturates quite rapidly, while that of the top remains slightly saturated and that of the bottom always dry, figure 14.

Drainage in a heterogeneous column (dry barriers)
This is an example of covers with capillary barrier effects from the authors' work in [26]. An initially saturated multilayer medium is drained. The layers are of materials with highly contrasted hydraulic properties. The values of the parameters characterizing the covering materials are:

Sand
Unreactive mining rejection The geometry is a vertical column of height = 110 cm containing a layer of unreactive mining rejection (fine silty material) 60 cm thick confined between two layers of sand (coarse material) 30 cm thick for the bottom layer and of 20 cm thick for the top layer. The medium was initially saturated (the interstitial pressure head was zero everywhere), thereafter a free drainage condition was imposed at the base of the column. The simulation is 60 days. Both layers of sand drain very quickly, figure 15. On the second day, the water content and pressure curves stabilized. The volumetric water content of the fine layer remains substantially at its saturation value, figure 16.

Conclusion
A finite element method has been developed for calculating the distribution of the interstitial pressure head and the volumetric water content in a water flow through a heterogeneous porous medium such as soil with variable saturation and initially dry. In the resolution, boundary conditions close to the actual conditions are taken into account. In particular, mixed-type conditions corresponding to atmospheric conditions on the soil surface. Conventional flow tests in the columns have proved the robustness and stability of the method with good conservation of the mass.