Contribution to Numerical Modeling of Water Flow in Variably Saturated, Heterogeneous Porous Media Contribution à la résolution numérique d’écoulement en milieux poreux hétérogènes variablement saturés

Abstract

Solutions to the Richards equation for water flow in variably saturated porous media are the focus of this paper. Working with field conditions, the extreme variability and complexity of soil, initial and boundary conditions can make the flow problem difficult to solve. This paper proposes to improve the computational efficiency of the mixed hybrid finite element (MHFE) method coupled with the variables transformation. The transform variables were introduced in order to simulate problems with convergence difficulty attributed to the presence of sharp wetting fronts. Furthermore, for better convergence behaviour, a technique that switches between the mixed-form and the pressure-head based form of the Richard’s equation was applied. Special attention was given to the top boundary conditions dealing with ponding or evaporation problems. In order to avoid non-physical oscillation problems, a mass condensation scheme was implemented in the model. Performance indicators in time and error of different options of the numerical model are defined, analyzed and classified. Thus, for each test case, a suitable numerical method that identifies which form of the Richards equation is best suited, the relevance of the switching technique as well as the utility of the transformation of the primary variable is possible. The results for the 1D Numerical test cases that have been performed matched those from the literature results. For evaporation and infiltration problem’s, the number of iterations needed to get the solution decrease when using the method of transformed pressure. Finally, knowing the soil heterogeneity, initial and boundary conditions, an agglomerative hierarchical clustering allows to analyze the need or not to transform variables and to use other options.

Key words:

• Modeling,
• numerical method,
• porous media,
• Richards equation,
• unsaturated,
• variables transformation

Résumé

Cet article est une contribution à la résolution numérique des écoulements en milieux poreux variablement saturés. La grande variabilité et la complexité des caractéristiques des sols, des conditions initiales et des conditions aux limites rendent les problèmes d’écoulements en milieu poreux difficiles à résoudre dans les limites acceptables de précision et de temps de calcul numérique. De nombreux efforts ont été consacrés dans la littérature récente sur le développement de solutions numériques. Ce papier propose d’améliorer l’efficacité de calcul de la méthode des éléments finis mixtes hybrides couplée à la transformation des variables. La transformation de variables est introduite dans le but de résoudre les difficultés de convergence dues à l’avancée d’un front raide de saturation en eau au sein du milieu poreux. Une attention particulière est portée aux conditions de saturation en eau proche de la surface du sol, telle la présence de flaques d’eau ou les phénomènes d’évaporation. Dans le but d’éviter des oscillations non physiques, un schéma de condensation de la masse a été implémenté, ainsi qu’une technique de basculement permettant de passer de la forme mixte de l’équation de Richards à la forme en pression uniquement. Des indicateurs de performances en temps et précision de calcul sont définis, analysés et classifiés. En conséquence, les options numériques optimales sont identifiées pour chaque cas testé. Par ailleurs, les résultats obtenus à l’aide des cas tests monodimensionnels tiennent très bien la comparaison avec ceux de la littérature. Connaissant l’hétérogénéité d’un sol, les conditions initiales et aux limites du domaine d’étude, une analyse par partitionnement des données propose de déterminer dans quels cas la transformée des variables ou les techniques de basculement entre les différentes formes de l’Equation de Richards sont adaptées.

Mots clés:

• Équation de Richards,
• modélisation,
• méthodes numériques,
• milieux poreux,
• non saturés,
• transformation de variables

1. Introduction

The unsaturated zone is by definition a multiphase system with at least two fluid phases present: air and water. An important assumption in deriving the Richards equation (Richards, 1931) is to assume that the air present in the unsaturated zone has infinite mobility. In other words, the air-phase pressure is assumed constant and equal to the atmospheric pressure, and air moves without interfering with water and/or contaminant. This assumption is reasonable in most cases because the mobility of the air phase is much larger than that of the water, due to the viscosity difference of the two fluids. Other assumptions in deriving Richards equation are constant water density and negligible porosity changes.

Taking into account these assumptions, the theory for transient flow of water into unsaturated-saturated soil is well understood and has been developed to a large extent in terms of solution of the nonlinear Richards equation. Vanderborght (2005) gives an overview of analytical solutions that can be found, for simple initial and boundary conditions, to define benchmark scenarios to check the accuracy of numerical solutions of the flow and transport equations. Working with field conditions, the extreme variability and complexity of soil characteristics, initial conditions, and varying boundary conditions can make the flow problem difficult to solve within acceptable limits of accuracy and computational effort. In view of this, most efforts have in recent years been concentrated on seeking numerical solutions (Fahset al., 2009; Mousaviet al., 2011; Serrano, 2004; Tian and Hu, 2007). Some studies to establish a comparison between numerical methods including the finite difference method, the finite element method, the finite volume method and the mixed finite element method have been published (Belfort and Lehmann, 2005; Farthinget al., 2003; Rees et al., 2004).

The term mixed method is used in problems where two or more physical variables are involved: in civil engineering both stress and displacement fields are approximated as primary variables. For fluid dynamics, the mixed finite element method allows us to compute an approximation of pressure head and velocity field simultaneously, e.g. with the same order of convergence. In order to obtain a positive definite matrix, hybridization is applied. In this paper, the solution has been approximated using the mixed hybrid finite element (MHFE) method of discretization (Arnold and Brezzi, 1985; Chavent and Jaffre, 1986). This technique is particularly well adapted to the simulation of heterogeneous flow field (Mosé et al., 1994; Nayagum, 2001; Youneset al., 1999). The reader interested by MHFE application for solving transport problems could refer to Wankoet al. (2009; 2011).

The objective of this paper is to improve the computational efficiency of the mixed hybrid finite element formulation for solving variably saturated 2D flow problems. Improvements were made by transforming the primary variables and applying a technique that switches between the mixed-form and the pressure-head form of the Richards equation. Special attention was given to the top boundary conditions dealing with ponding or evaporation problems. A mass condensation scheme was used in order to avoid oscillation problems related to the discrete expression. Finally, using statistical analysis on time-indicator and error-indicator parameters, suitable numerical options were proposed for different types of flow problems including soil heterogeneity, different initial and boundary conditions.

2. Materials and methods

2.1 Theory

2.1.1 An overview of model approaches for water flow in variably saturated porous media

According to Köhneet al. (2009), classical physically based model concepts for water flow in structured soils can be broadly classified into three categories.

2.1.1.1 Continuum model

The traditional conceptualization of soil is that of a porous medium with continuous properties. Variably saturated water flow through such a porous system is uniform and at local equilibrium. The simplest and most widely applied mechanistic, Darcy-scale model for water flow in the unsaturated zone consists of the Richards (1931) equation. This model was used here.

2.1.1.2 Bi- or multi-continuum model

The mobile–immobile model assumes that the soil micropore network is poorly connected and of such a low permeability that water is immobile in the vertical direction. Whereas the dual-permeability model assumes that the porous medium consists of two overlapping pore domains, with water flowing relatively fast in one domain and slow in the other domain. In practice, this approaches assume a porous system which is permanently in the shrinkage state (Coppolaet al., 2012). The continuum models cannot explain the physical processes at pore scale.

2.1.1.3 Network model

This model allows quantification of flow when a spatial averaging approach based on a representative elementary volume cannot be considered in the study area. As an example, the discrete fracture models is well adapted to modelling flow and transport in fracture network (BERKOWITZ, 2002). The pore network model is another main concept for this model (JOEKAR-NIASAR et al., 2012).

2.1.2 The Richards equation

The below formulation (1) physically describes the flow in a variably saturated porous medium by taking into account the elastic storage effects.

(1)

Where:

• f(x,z,t) is the sink/source terms [T^‑1];

• x and z (depth) are the spatial coordinates [L];

• t is time [T];

• C(h) is the soil moisture capacity [L^‑1];

• K is the unsaturated hydraulic conductivity [LT^‑1];

• h is the soil water pressure head [L];

• S_s is the specific storage [L^‑1];

• S_w is the degree of saturation (-).

Despite the development of fairly versatile numerical simulators for the Richards equation, there is still a deep interest in the development of exact and approximate closed-form solutions of the Richards equation, particularly for flows in heterogeneous media, which is the main essential characteristic of real soil. Numerical modeling of the Richards equation has yielded significant progress, as has the development of faster, cheaper computers. Some very useful results have been obtained. As stated above in the introduction, MHFE method has shown its effectiveness in solving the Richards equation.

2.1.2 Mathematical formulation and MHFE implementation

A two-dimensional domain Ω is defined and space-discretized into triangular elements G. The continuity of pressures and fluxes between two adjacent elements is valid for all the interior edges E_i (∀i = 1,2,3) of the domain Ω. The average water pressure by edge is chosen as the unknown of the symmetric positive definite linear system to solve. The Darcy flux is approximated over each element by a vector belonging to the lowest order Raviart-Thomas space (RAVIART and THOMAS, 1977). On each element the vector function has the following properties: is constant over the element G, is constant over the edge E_i of the triangle, ∀i = 1,2,3, where is the normal unit vector exterior to the edge E_i. is perfectly determined by knowing the flux through the edges (CHAVENT and ROBERTS, 1989).

(2)

Where Q_G,Ej is the water flux over the edge E_j belonging to the element G [L² T^‑1] and the vector fields basis used as basis functions over each element G. These vector fields are defined by , where δ_i,j is the Kronecker symbol:

2.1.3 Variables transformation

PAN and WIERENGA (1995) presented a method to simulate problems with convergence difficulty attributed to the presence of sharp wetting fronts. The pressure head variable h is then transformed into new dependent variables ĥ (3a) and K̂ (3b):

(3a)

Where:

• κ is a constant (-0.05 < κ < -0.01 cm^‑1) independent of both the K(h) and C(h) relationship (Pan and Wierenga, 1995);

• he is a free parameter, referred as the air entry value [L] (Ippisch et al., 2006);

• ĥ is the transformed pressure head [L].

The transformed hydraulic conductivity is defined as K̂ :

(3b)

Where:

• K̂ is the transformed hydraulic conductivity [LT^‑1].

2.1.4 The transformed mass balance equations

Based on the approximation of variables over each element of discretization G, the standard pressure based form of the Richards equation (HILLEL, 1980) can lead to large mass balance errors (4a); while the mixed form (RICHARDS, 1931) has improved properties with respect to accurate mass conservative solutions (4b), but it can have convergence difficulties for dry initial conditions.

(4a)

∀G over a domain Ω, for t ⊂ ]0,T[

(4b)

∀G over a domain Ω, for t ⊂ ]0,T[

Where:

• is the substitution variable of C, the specific water moisture capacity [L^‑1];

• is the substitution variable of θ, the volumetric water content.

The average hydraulic conductivity for each element is estimated using the relationship K_G = Ks_GKr_GK^A_G. Where the sub index G denotes the element, Ks is the saturated hydraulic conductivity [LT^‑1], K^A is a dimensionless anisotropy tensor, Kr is the relative hydraulic conductivity function (dimensionless). Kr_G is considered to be the maximum value among the three Kr_E values, which are calculated using the values of pressure at the interior edges and the modified Mualem-van Genuchten expression (IPPISCH et al., 2006).

2.1.5 Switching technique

The primary variable switching technique has proved to be an effective solution strategy for unsaturated flow problems (DIERSCH AND PERROCHET, 1999; HAO et al., 2005). It is unconditionally mass conservative. Better convergence behaviour is achieved using this technique, compared to both the mixed-form and pressure-head based form of the Richards equation. The primary variable is switched at each iteration, using the following criterion: if then pressure head is used as primary variable, if not then a mixed-form of the Richards equation is used. The tolerance for this switching procedure, tol_f, is predefined by the user (0 ≤ tol_f, ≤ 1).

2.1.6 Mass condensation scheme (Mass lumping)

BELFORT (2006) proposed the use of a mass condensation scheme in order to avoid unphysical oscillations, which are related to the time-dependant terms appearing in the diagonal coefficients of the mass matrix used for the numerical solution. An expression of the flux at each edge Q_G,Ei is defined in terms of steady state and transient flow regimes, for details see (YOUNES et al., 2006).

2.1.7 Top boundary conditions

VAN DAM and FEDDES (2000) developed a procedure for 1D model that gives special attention to the top boundary condition, which is important for simulations with ponded water layers or with fluctuating water levels close to the soil surface. This procedure (Figure 1) switches from head to flux controlled boundary condition and vice versa. This algorithm for the management of boundary conditions was adapted to be used with the mixed hybrid formulation and 2D flow problem.

Figure 1

Procedure to select head or flux top boundary condition. Adapted from Van Dam and Feddes (2000)

Procédure de sélection de la charge ou du flux imposé à la limite supérieure. Adapté de Van Dam and Feddes (2000)

2.1.8 Mass balance error

Following the approach of CELIA et al. (1990), a mass balance measure MB (5) is defined in order to test the ability of the model to conserve mass. The accuracy of the numerical scheme is evaluated by computing a global mass balance error ε_MB (6):

(5)

(6)

Where:

• Total additional mass in the domain =

• Total net flux into the domain =

• nt is the total number of time steps and

• nm is the total number of elements in the domain.

Additionally, a mass balance ratio can be computed for each time step MBⁿ (7), with their corresponding error ε_MBⁿ = 1 − MBⁿ.

(7)

E_1G (8) and E_2G (9) are local mass balance error computed from edge and mesh properties values respectively.

(8)

With E_1Gi the local mass of edge i

(9)

A low global mass balance error is a necessary but not a sufficient condition to ensure accuracy of the solution. Mass balance error can decrease even if the solutions do not converge (KOSUGI 2008; TOCCI et al., 1997). KOSUGI (2008) concluded that is important to check the mass balance and solution convergence at each time step when using the discretization scheme of CELIA et al. (1990) to simulate unsaturated water flow.

2.1.9 Maximal convergence errors for pressure head and water content

Maximal convergence errors for pressure Δh_max (10) and water content Δθ_max (11) are calculated at the end of each time step as follows:

(10)

(11)

where max is the generator of the maximal value in the entire spatial domain, m is the iteration index and n is the time step index.

2.1.10 Discrepancy between average pressure and arithmetic mean of edge pressures

Δh_dmax (12) represents a measure of disagreement or discrepancy between the average pressure calculated in the element using the mixed hybrid formulation and the arithmetic mean of the three edge pressures belonging to this element.

(12)

where max is the generator of the maximal value in the entire spatial domain and Th_{G, Ei} are the three water pressure traces of the element G.

2.1.11 Numerical solution and convergence criteria

The numerical method will lead to a system of linear equations, where the unknowns are the water pressure traces  (Th). The number of unknowns is equal to the number of edges to which the pressure has not been imposed. The matrix associated with the hydrodynamics equations system is symmetric and definite positive. Therefore, it can be effectively solved by the conjugate gradient method, preconditioned with an incomplete Cholesky decomposition using the Eisenstat procedure (Eisenstat, 1981). The iteration process for the hydrodynamic calculation when using the mass condensation scheme is stopped when the following conditions are met:

• The difference between the calculated values of edge pressure head between two successive iteration levels is smaller than an absolute iteration convergence tolerance predetermined by the user (SHAHRAIYNI and ASHTIANI, 2009).

• The iteration convergence test, which involves both absolut (equations 10 and 11) and relative error for pressure and water content is satisfied (KAVETSKI et al., 2001).

• The difference between the calculated values of water content between two successive iteration levels is smaller than a tolerance predetermined by the user (HUANG et al., 1996).

2.1.12 Time Control

The temporal space ]0,Ts[ is discretized in temporal increments (Δt), which are automatically adjusted at each time level according to the following rules:

1. There is a minimum and a maximum time step (Δtmin and Δtmax) that are specified by the user. So then, Δtmin ≤ Δt ≤ Δtmax.

2. A maximum temporal increment allowed is estimated by setting as a maximum value 0.5 for the dimensionless Fourier number (Fn), represented by a ratio of the Courant-Friedrichs-Lewy (Co) and Peclet (Pe) numbers by using the Kirchhoff transformation (EL-KADI and LING, 1993). The maximum allowed Δt for the hydrodynamics (Δmax_{hydrodynamics}) will be the minimum value in the entire spatial domain.

3. The initial time step Δt will be equal to the minimum value between Δtmax_hydrodynamic and a smaller defaut initial value (Δt_init). For the next time levels, a heuristic method (BELFORT, 2006; SIMUNEK et al., 2005) will be used but always respecting the rules 1 and 2.

2.2 Numerical test cases

2.2.1 Pan and Wierenga (1995) flow problem

Twelve one-dimensional cases (PAN and WIERENGA, 1995) were simulated using layered or uniform soil profiles (Figure 2). Soil properties are listed in table 1. The initial pressure and boundary conditions applied in each case are described in table 2.

Figure 2

Soil profiles for Pan and Wierenga [1995] test cases. Layered soil profile (left), uniform soil profile (right)

Profils de sols issus des cas tests de Pan et Wierenga [1995]. Sol multi-couches (gauche), sol uniforme (droite)

Table 1

Soil parameters used in the modified Mualem and van Genuchten Model

Paramètres de sols d’après le modèle modifié de Mualem and van Genuchten

Table 2

Initial pressure and boundary conditions

Conditions initiales et aux limites

A transformed pressure (PAN and WIERENGA, 1995) was introduced as the dependent variable with the MHFEM and results were compared to those without using transformation of variable. Moreover calculations were performed using the Richards equation on their h-based form, mixed-form or using a switching method between these two forms. Mixing these entire variables, a total of 72 simulations have been performed.

2.2.2 Ponded conditions at a dry soil and high evaporation at a wet soil

The ability of the model to deal with problems related with the top boundary condition was tested for two cases of extremes conditions at sand soil. Soil characteristics are given by soil 4 in table 1.

For ponded conditions, the rainfall rate was 1 000 mm∙d^‑1. The initial conditions on water content were equal to 0.1. At the reference case, the hydraulic head gradient at the soil surface is large enough to absorb the infiltration rate of 1000 mm∙d^‑1, until a time of simulation t = 0.008 d. The cumulative amount of infiltration obtained in the reference case was 39 mm.

Concerning the other test case, the potential evaporation rate was 5 mm∙d^‑1. The initial pressure was -2000 mm and the cumulative actual evaporation obtained in a period of five days at the reference case was 11 mm.

2.3 Agglomerative Hierarchical Clustering (AHC)

The Agglomerative Hierarchical Clustering (AHC) procedure is used to make up homogeneous groups of objects (classes) on the basis of their description by a set of variables.

AHC was applied here to classify in homogeneous groups on the one hand all the indicator parameters included in table 3, on the other hand the different numerical options. The statistical analysis was performed with the software XLSTAT 2010.

Table 3

Computational parameters

Paramètres de calcul numérique

2.3.1 AHC Variables definition

For this statistical analysis, two types of variables were defined:

2.3.1.1 Discrete variables

According to the following definitions (Table 3), discrete variables describing the choice of the model is composed by five alphanumeric characters identifying: the type of boundary condition (N or D), the specified top boundary condition (1, 2, 3 or 4), the initial condition on pressure (A, I, or W), the tolerance used for the switching procedure representing the type of the Richards equation used (H, M or S), and the type of primary variable used (T or P).

For example: N1AHT means the simulation results are obtained by setting a Neumann boundary condition (N), with a high flux imposed (1), soil is (A). The Richards equation was solved using the h-based form (H) and the primary variable used was a transformed pressure head (T). There are in total 72 discrete variables.

2.3.1.2 Continuous variables

Continuous variables carry the quantitative information obtained after simulation. There are 15 continuous variables (time and error indicators). See table 4 on computational and description parameters.

Table 4

Definition of discrete variables

Definition des variables discrètes

3. Results and discussion

3.1 Pressure head and water content distributions

Results show that the method is numerically robust for all cases of variably saturated, heterogeneous media, and first or second type boundary conditions. Pressure head and water content distributions were in very good agreement with PAN and WIERENGA (1995) results (Figure 3).

Figure 3

Pressure head and water content distributions

Distributions de la charge de pression et de la teneur en eau

Figure 3 (continuation)

3.2 Suitable method for hydrodynamic simulation in porous media

3.2.1 Statistical analysis

Table 5 summarizes the general statistical analysis of the data:

Table 5

Statistical analysis of the data

Analyse statistique des données

3.2.2 Clustering results

Data analysis enabled the re-grouping of discrete variables, or continuous variables, in homogeneous groups.

3.2.2.1 Re-grouping of discrete variables

In order to re-group the discrete variables, figure 4 shows the dendrogram representing the hierarchy obtained using a euclidean distance as the dissimilarity metric between points and the Ward Agglomeration method. Three class centroids were distinguished (Table 6 and Figure 4). It appears that class  2 has the minimum continuous variables NI, E_1G max, E_2G max, ∆h max, ∆θ max and ∆h_d max (Table 6).

Figure 4

Dendrogram for discrete variables

Dendrogramme des variables discrètes

Table 6

Class centroids

Centre des classes

The central objects for class 1, class 2 and class 3 are represented by discrete variables D3AHT, N1WST and D3WSP respectively. So that it is possible to make link between the choice of the numerical option for the model (discrete variables) and their advantageous characteristics. Distances between the centroids of class 1 to class 2 and class 3 are in the same order of magnitude (9.9 x 10⁴). The distance between the centroids of class 2 and class 3 is shorter (5.3 x 10³).

From the results above it can be deduced that:

• The first class is distinguished by arid soil as initial conditions. All simulations concerning this type of initial pressure condition are in this group, regardless of the boundary conditions that were imposed, the form of the equation to solve, or the primary variable used (Figure 4).

• The third class is characterized by Dirichlet boundary conditions with a positive pressure imposed, except for soil arid initial conditions (Figure 4).

• The second class deals on the one hand with Neumann boundary conditions except for arid soil initial conditions, and on the other hand with Dirichlet boundary conditions with negative pressure imposed, except for arid soils initial conditions (Figure 4).

• The three classes are regardless of the form of the equation to solve and the primary variable used.

3.2.2.2 Re-grouping of continuous variables

In order to re-group the continuous variables, figure 5 shows the dendrogram representing the hierarchy obtained using an euclidean distance as the dissimilarity metric between points and the Ward Agglomeration method. Three class centroids were distinguished. The central objects for class 1, class 2 and class 3 are presented by discrete variables Time steps, ∆t min, and E_1G max respectively. Distances between the centroids of class 2 to class 1 and class 3 are in the same order of magnitude (2.5 x 10⁵). The distance between the centroids of class 1 and class 3 is shorter (3.7 x 10⁴).

Figure 5

Dendrogram for continuous variables

Dendrogramme des variables continues

From the results by class it can be deduced that:

• The third class represents only global and local mass balance errors. Therefore, it can be considered as an indicator of mass conservation.

• Pe max and Co max numbers are numerical parameters involved in the time stepping size procedure. As it was expected, they are directly related with explicit time parameters (∆t_av, ∆t_min, ∆t_max) and they were grouped together in class 2.

• The parameters ∆h_max, ∆θ_max and ∆h_dmax are indicators of precision, which are directly related with the total number of iterations and time steps. They were grouped in class 1.

The previous comments associated with the fact that the shortest distance between centroids was found between class 1 and class 3, suggest that the global and the local mass balance errors have more proximity to the precision parameters (class 3) than to the time controlling parameters (class 2).

3.2.3 Selection of suitable numerical options

The re-grouping of the qualitative variables by classes enables the definition of centers of gravity for each class. In the following paragraphs we propose which simulations are the best suited for each test case. There were six numerical options for each test case.

The choice of the models was performed by sorting in ascending order with a progressive constraint the observation values of the three quantitative variables constituting the centers of gravity of each class. That is to say, a number of time steps minimum, a ∆t min maximum and the smallest E_1G max. Thus, for each test case, a suitable numerical method that identifies which form of the Richards equation is best suited, the relevance of the switching technique as well as the utility of the transformation of the primary variable is possible.

The test cases 1.1 and 2.1 simulate the infiltration in an arid soil with a Neumann condition, by imposition of low and high flux, respectively. The best adapted model proposes the mixed form of the Richards equation with transformation of the variable of pressure. It should be noted in both cases, that the first three most relevant options used the variables transformation technique.

The test cases 1.2 and 2.2 simulate the infiltration in an intermediate soil with a Neumann condition, by imposing low and high flux, respectively. The best adapted model proposes the application of the variables transformation technique for both cases and the application of the switching technique for test case 1.2 and h-based form of the Richards equation for test case 2.2.

The test cases 1.3 and 2.3 simulate the infiltration in a wet soil with a Neumann condition, by imposing low and high flux, respectively. The best adapted model proposes for test case 1.3 to not transform the pressure variable and to apply the switching technique. For test case 2.3, the transformation of variables and the h-based form of the Richards equation were proposed.

The test cases 3.1 and 4.1 simulate the infiltration in an arid soil with a Dirichlet boundary condition, by imposing positive and negative pressures, respectively. The best adapted model proposes the transformation of the variable of pressure coupled to the switching technique for test case 3.1 and the non transformation of pressure coupled to the mixed-form of the Richards equation for test case 4.1.

The test cases 3.2 and 4.2 simulate the infiltration in an intermediate soil with a Dirichlet boundary condition, by imposing positive and negative pressures, respectively. The best adapted model proposes the mixed-form of the Richards equation for both cases, the transformation of pressure for test case 3.2, and the non-transformation of the primary variable for test case 4.2.

The test cases 3.3 and 4.3 simulate the infiltration in a wet soil with a Dirichlet boundary condition, by imposing positive and negative pressures, respectively. The best adapted model proposes the mixed-form of the Richards equation and the transformation of the primary variable for both test cases.

From this analysis, it is deduced that the less indicated models, according to the established criteria of selection, are those applying the non-transformation of the primary variables coupled to the mixed-form or the h-based form of the Richards equation, with two exceptions (N1WMT, D4WHT). In particular, the models coupling the h-based form of the Richards equation and the non-transformation of the primary variable are the less indicated for the problems applying Neumann boundary conditions with a low flux imposed. H-based form of the Richards equation is less indicated for problems applying Dirichlet boundary conditions with a negative pressure imposed, while the non-transformation of the primary variable are the less indicated for problems applying Dirichlet boundary conditions with a positive pressure imposed .

3.3 Extreme top boundary conditions

Table 7 shows the different values calculated as cumulative amount of infiltration and the time step sizes for each simulation while the cumulative actual evaporation obtained in a period of five days are presented in table 8.

Table 7

Total infiltration

Infiltration totale

Table 8

Total evaporation

Évaporation totale

The number of iterations needed to get the solution decrease when using the method of transformed pressure. Cumulative amount of infiltration is nearly the same for all the methods (Figure 6). The selection of an appropriate equivalent conductivity when simulating infiltration in dry soil or high evaporation from wet soils is important (Figure 6). Geometric, weighted and integrated formulations produce better solutions than a traditional scheme using a mean conductivity calculated with a mean pressure head (BELFORT and LEHMAN, 2005). However, the use of the geometric mean to estimate the hydraulic conductivity underestimates the water fluxes or leads to convergence problems. The method proposed here to estimate the mean hydraulic conductivity for an element consists in assigning the maximum value calculated with the edge pressure heads. As it can be seen, results show good agreement with those presented by VAN DAM and FEDDES (2000).

Figure 6

Illustration of infiltration rate (left) and evaporation rate (right)

Illustration des taux d’infiltration (gauche) et des taux d’évaporation (droite)

Conclusions and Outlooks

Mixed hybrid finite element was used to simulate different flow problems in variably saturated porous medium. This formulation is based on Raviart-Thomas’ space properties. Different techniques were used, such as the transformation of the primary variable of pressure. For a better convergence behavior, a technique that switches between the mixed-form and the pressure-head based form of the Richard’s equation was applied. Special attention was given to the top boundary conditions dealing with ponding or evaporation problems. In order to avoid non-physical oscillation problems, a mass condensation scheme was implemented in the model.

The model was checked by performing several simulations and the comparison of the numerical approximations obtained with reference results found in the literature. A set of test cases were simulated by using different model options, depending upon the primary variable used, the form of the equation to solve, and the initial and boundary applied conditions. As simulation result, a group of time and error indicator parameters was defined. Observations of these parameters were analyzed statistically, and re-grouped according to their dissimilarities. Their re-grouping according to dissimilarities gave as result information that could be very useful for the selection of the most suitable option-model to apply, depending on the initial and boundary conditions to simulate. Main results are summarized below:

• The mixed form of the Richards equation with transformation of the pressure variables are suitable both for arid soil with Neumann condition and infiltration in a wet soil with Dirichlet condition;

• The h-based form of the Richards equation with transformation of the pressure variables are appropriate for the high flux infiltration in a wet soil;

• The transformation of the pressure variables coupled to the switching technique suitable for infiltration in a arid soil with positive Dirichlet condition;

• Simulation of infiltration in a very dry soil or high evaporation in a wet soil is made possible by a suitable choice of equivalent conductivity.

Hence modeling infiltration in variably saturated porous media with different type of boundary conditions could be done with appropriate numerical options.