Mathematical Problems in EngineeringVolume 2012 (2012), Article ID 187634, 23 pagesdoi:10.1155/2012/187634

Research Article

The Finite Volume Formulation for 2D Second-Order Elliptic Problems with Discontinuous Diffusion/Dispersion Coefficients

Stefano Ferraris, Ivan Bevilacqua, Davide Canone, Davide Pognant, and Maurizio Previati

Dipartimento di Economia e Ingegneria Agraria, Forestale e Ambientale (DEIAFA) Sez. Idraulica, Via Leonardo da Vinci 44, 10095 Grugliasco , Italy

Received 18 July 2011; Revised 2 November 2011; Accepted 17 November 2011

Academic Editor: Moran Wang

Copyright © 2012 Stefano Ferraris et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We propose a finite volume method for the numerical resolution of two-dimensional steady diffusion problems with possibly discontinuous coefficients on unstructured polygonal meshes. Our numerical method is cellcentered, secondorder accurate on smooth solutions and based on a special numerical treatment of the diffusion/dispersion coefficients that makes its application possible also when such coefficients are discontinuous. Numerical experiments confirm the convergence of the numerical approximation and show a good behavior on a set of benchmark problems in two space dimensions.

1. Introduction

In the last decades, finite volume methods have been greatly successful in solving engineering models of flows in porous media on complex geometrical domain because the finite volume formulation works on general polygonal and polyhedral meshes. This great mesh flexibility is now combined with a strong theoretical foundation, that is, convergence analysis and error estimates are available on the simplest mathematical models [1]. In this work, we address a version of the finite volume method that is popularly known as “the diamond scheme” and was originally presented for the advection-diffusion equation in two dimensions in [2–4] and successively extended to convection-dominated problems in [5, 6] and nonlinear flow problems in partially saturated porous media [7].

The numerical treatment of both diffusion/dispersion flux and of the partially saturated flow is based on averaging one-sided gradients and diffusion (conductivity) tensors. Many papers addressed the issue of conductivity averaging both in finite differences and in finite elements. A seminal work [8] started with a finite difference application to the unsaturated flow equation. Later on they, dealt with the problem in a finite volume scheme [9]. Recent reviews on this topic can be found in [10–12]. They also made performance comparisons using different averages. The latter paper shows in conclusion that their scheme, namely, an averaging schemes based on Darcian mean principle used in the framework of either vertex-centered or cell-centered approach compare favorably to other methods for a range of test cases. They compared their method with the one of [11]. Aside of the numerical comparisons, in [13] it is remarked that the usage of the harmonic average can be motivated on the basis of physical arguments.

Another reason for dealing with the tensor averaging issue is related to the scale problems. Typically, the Darcy equation is written at a scale much smaller than the scale of practical interest [14, 15]. This gives boost to new numerical techniques to able to deal with the heterogeneity of hydrodynamic parameters. One possibility is to obtain either analytically or empirically explicit equations for the scale of interest, eliminating other scales in the problem [16].

Several methods have been recently developed, such as, the heterogeneous multiscale method (HMM) [16] and the multiscale finite volume element method [17]. Babuška in the 70s motivated the multiscale finite element method [18–20]. Multiscale methods have been proposed also in the case of saturated flows in heterogeneous porous media and also applied the multiscale finite volume element method to the Richards equation [21, 22]. In [14], a multiscale method based on the framework of HMM was recently extended to solve the Richards equation with random coefficients.

Another feature of their numerical scheme is the formula estimating the macroscopic flux in which the unsaturated hydraulic conductivity can be calculated as a diagonal tensor. In particular, it is worth noting that a finite volume high-order accurate approximation of the pressure head field also allows one to achieve a better resolution in the approximation of other important fields like the components of the hydraulic conductivity tensor at the mesh vertices [7]. For the sake of simplicity, in this work it was considered the steady flow in a layered porous medium in the presence of a source term. It is the Poisson equation which is suitable to test new schemes which will be applied in further work to the more general unsaturated transient cases.

We propose a novel technique that automatically adjusts when the diffusion tensor is discontinuous across a mesh interface shared by two adjacent cells. This technique is general and can be easily implemented in any finite volume scheme that has an explicit numerical flux and may result in a particularly efficient strategy for the numerical resolution of problems where both the diffusive/dispersive and convective phenomena are simultaneously significant. In fact, in such problems high-resolution finite volume methods are normally coupled with mixed finite elements following the criterion of choosing the best available technique in accordance with the nature of the equation to be discretized. In fact, the method is more suitable to the numerical discretization of the diffusive part of the model, and the finite volume method gives an accurate and stable discretization of the convective part of the model, even in presence of strong convective fields. This approach was proved successful in the numerical modeling of oil reservoir problems [23] and of groundwater flow and transport of contaminants in porous media, combared with [24–26]. Moreover, other different engineering areas may benefit from this new technique, such as [27–38].

The outline of the paper is as follows. In Section 2, we introduce the mathematical model and discuss the formulation of the finite volume methods based on vertex reconstruction. In Section 3, we present the numerical treatment of the diffusion tensor. In Section 4, we confirm the theoretical results with numerical experiments. In Section 5 we offer final remarks and conclusions.

2. The Mathematical Model and the Finite Volume Formulation

Let be a polygonal domain with boundary .

We consider the steady diffusion problem for the scalar solution field given by where is the diffusion tensor describing the material properties, is the forcing term and the boundary function that defines nonhomogeneous Dirichlet conditions on the boundary . Under suitable regularity assumptions on , , , and , it turns out that the diffusion problem is mathematically wellposed, a unique solution exists [39] and such a solution is continuously dependent on the model data. In particular, due to the elliptic nature of the model, the diffusion tensor is normally represented by a strictly positive definite matrix for every in . The components of the diffusion tensor may be discontinuous in the computational domain; in such a case and without loss of generality, we assume that the mesh is conforming with each discontinuity of .

The numerical approximation to (2.1)–(2.2) is performed on a sequence of polygonal partitions of the domain . For any mesh , the subscripted label is the mesh size and is defined by

where is the diameter of the polygonal cell . We also label the numerical solution calculated on the mesh by . Also,(i)a generic mesh vertex is denoted by and its coordinate vector by . We will also find it convenient to introduce a local numbering of the vertices, for example, , and to ease notation to denote the vertex position of the -th vertex by (instead of );(ii) a generic cell interface or a boundary edge is denoted by , its center (i.e., its edge midpoint) by , and its measure (edge length) by ;(iii) a generic polygonal cell is denoted by its measure (area) by , its center of gravity by , and its boundary by .

The orientation of each mesh interface is reflected by its unit normal vector , which is fixed once and for all. For any mesh cell and any face of the polygonal boundary , we define the unit normal vector that points out of and we also use the notation .

Let the algebraic vector be the numerical solution, where each approximates the cell average of the scalar solution over the cell . The finite volume scheme for on the computational mesh reads where is an evaluation of at face taken from the side of ; is the discrete face gradient built inside the diamond cell centered at face ; is the geometric vector perpendicular to , pointing from to , and with lenght ; is the cell average of the right-hand side term over the cell :

Let and be first-order approximations of the diffusion tensor within the control volumes and . For example, we can take and , where and are internal points (not necessarily the barycenters) of and , respectively. Now, we can consider either inside and inside or a suitable mean of these two tensors. In particular, we can deal with the arithmetic mean:

or with the harmonic mean:

Both approaches deserve a special care when is discontinuous across and a proper definition of the discrete gradient is also needed in order to preserve the property of flux conservation.

To define the numerical gradient, a special control volume is built around this interface, which has a quadrilateral shape in two dimensions and is named “diamond cell”. The geometry of a diamond cell is shown in Figure 1, which plots the diamond cell centered at the mesh face with vertices and and shared by the cells with centers and . The diamond cell can also be seen as the union of the subdiamonds and , which are the triangular cells sharing as common base and having, respectively, vertices and , the centers of gravity of cells and . The four vectors , , , and shown in Figure 1 are respectively orthogonal to the four boundary faces , , , and , and their lenght is equal to the lenght of the corresponding face. Instead, when is a boundary face, thus defined by ( being the boundary of the computational domain ), the diamond cell associated to coincides with , that is, it is the triangle defined by and the vertex .

Figure 1: Geometry of the 2D diamond cell.

The numerical diffusive flux is built by using a constant approximation of the solution gradient and an evaluation of the diffusion tensor within the diamond cell, as discussed before. Let us give the formulas for an internal mesh face , that we suppose to be shared by two cells and . All these formulas can be easily adapted to the case of a boundary face by taking equal to the center of gravity of face .

Using the Green-Gauss theorem yields: where is the unit vector orthogonal to and pointing out of . If the restriction of to the face of is an affine function, the boundary integral on in (2.8) only depends on the values of at the vertexes of and the constant vector provided by the formulas that we derive below must coincide with the gradient of .

The Gauss-Green theorem applied to the solution gradient integrated on the 2D diamond cell of Figure 1 yields: Let us introduce the vectors and , which are, respectively, orthogonal to the edge connecting to and to , and whose lenght are equal to the length of these edges. Using the geometric relations and into (2.9) and rearranging the terms yields:

To derive the gradient formula, we replace the function values and with the cell unknowns and , respectively, and the function values and with the corresponding nodal unknowns and . We obtain which is a constant approximation of the solution gradient within the quadrilateral cell .

With any vertex of the mesh , we associate the reconstructed vertex value , where denotes the subset of the mesh cells which share the vertex . The interpolation weights are assumed to verify the consistency relations [3]: Alternative constructions of the interpolation weights are available in the literature, see, for instance, [40]. The main advantage offered by these alternative choices concerns the control of positivity of the weights, which turns out to be significant when we are aimed at developing a numerical method with a maximum or minimum principle. However, as this topic is out of the scope of the current work, we will not pursue it anymore.

The interpolation weights are obtained by solving the reconstruction problem that approximates the cell-averaged data set by the affine function: on the covolume and in a least square sense, compared to [3, 41]. The reconstructed value at vertex is now given by taking . The coefficients are the minimizers of the least squares functional: Imposing the zero gradient condition, that is, , yields a linear system for the coefficients , whose solution returns the interpolation weights. For completeness of exposition, we briefly review the derivation of the weight formula given in [41]. Let ; a straightforward calculation yields:

Let be the number of cells of ; we reformulate the linear system written above as For convenience, we consider a suitable local numbering of the cells of , for example, , and we introduce the two matrices:

and the column vector , which collects the solution averages to be imposed. Using such definitions, the linear system (2.16)–(2.18) becomessince matrix is, obviously, symmetric and positive definite and, thus, nonsingular. The coefficient that provides the vertex value is given by We require that such a coefficient be an average of the data in with coefficients , that is, . By comparison with (2.21), we have the final weight formula:

Conditions (2.12) can be checked by a direct calculations.

The functional in (2.14) can be suitably modified to take into account different kind of boundary conditions, for example, Neumann or Robin, compared to[41]. The values at the vertexes on the Dirichlet boundary are constrained to the boundary data, for instance for a homogeneous condition. We also mention that this choice of weights provides a very robust technique, compared to the study of locking effect due to the accuracy of the reconstruction [42].

3. Treatment of Discontinuous Diffusion Tensors

We discuss, here, how to treat the case of a diffusion tensor that is discontinuous across an interface shared by the cells and , that is, . Since we suppose that is discontinuous across and the normal flux of the exact solution is continuous across , the normal component of the solution gradient must be discontinuous. Consequently, a numerical approximation of the diffusive flux across like

cannot be consistent whenever is some kind of average between and and the constant vector approximates   over all the diamond cell . To tackle this problem, we consider two distinct approximations of the solution gradient within the subcells and , denoted by and , respectively, and impose directly the flux conservation. To derive an expression for such one-sided gradients, we introduce an additional unknown along face and we apply the Gauss-Green Theorem as for the derivation of in the previous section. Also, we recall that and and we use a similar notation for the normal vector from the other side of , so that and . The two one-sided gradient formulas are

Now, we search for a tensor , which is an average of the diffusion tensors and , and a gradient vector , which is an average of the gradient vectors and , that makes it possible to express directly the flux conservation as Since we expect that the normal component of the vectors and be the same, we impose directly this condition by requiring that a scalar coefficient exists such that (recall that ). We also assume that be a weighted average of and . Let us be given two nonnegative coefficients and such that and such that

To this purpose, a number of possible choices have been proposed and widely investigated in the literature for and . For example, a popular option is to take the volume of the two subdiamonds and normalized by the total volume of the diamond, for example, and . The diffusion tensor can be chosen as the simple average of and , for example, , or as the harmonic average, for example, . Much more attention must be paid when is discontinuous across as none of the previous choices may provide the correct flux information across the common interface of cells and . Our aim in the next developments is at determining a value for and of the average diffusion tensor in terms of and that ensures flux conservation (3.3) for any given pair of coefficients and . Our derivation is similar to that considered in [43], and as pointed out therein, the resulting matrix is the same obtained for other purposes in the field of upscaling of conductivity, either by means of the large average techniques [44] or using the homogenization theory in the case of layered materials [45].

To this purpose, we first derive an expression for . We multiply (3.4) by and use flux conservation to obtain

from which we have Likewise, we multiply (3.4) by and use flux conservation to obtain

from which we have We multiply (3.8) by and (3.10) by , we sum the resulting relations and use (3.5) to have

Solving this last equation for gives us the formula: Now, we derive an expression for . In view of (3.5) and (3.3), it holds that First, in (3.13), we substitute the expression for provided by (3.4), we collect the factor and we obtain: Second, in (3.13), we substitute the expression for provided by (3.4), we collect the factor and we obtain Third, we multiply (3.14) by , (3.15) by , we add the two resulting equations, we use again (3.5) and (3.6), and we obtain Finally, in (3.16) we substitute the scalar factor given by (3.12), we collect the vectors and and we obtain: Equation (3.17) suggests us to set where

4. Numerical Experiments

In this section, we present and discuss a set of numerical results to show the convergence behavior of the finite volume method when we compute the diffusive flux using the technique described in the previous sections. In all the test case that we present in this section, we compare the behavior of the new numerical treatment that we consider in this paper with the weighted average, which is the standard approach in the method [5, 41, 42]. In fact, from the formulation given in equation (3.18), it follows that the average diffusion tensor is equal to the weighted average plus a correction term. This correction term is specifically designed to take care of possible discontinuities in the diffusion coefficients.

The finite volume formulation based on the least squares reconstruction of vertex values leads to a linear system for the cell-average unknowns whose coefficients matrix is generally unsymmetric although displaying a symmetric nonzero pattern. The positive definiteness of such system is still an open issue and this fact is also the major difficulty for the development of a full theory of convergence of such method. Theoretical results that prove coercivity are available only for meshes of slightly deformed quadrilaterals [2–4].

We solve such linear system using the MA41 routine of the HSL software collection, which implements an unsymmetric multifrontal sparse LU factorization technique especially designed for matrices with a symmetric nonzero pattern and unsymmetric values. Different software packages like UMFPACK and standard Krylov methods like BiCG-Stab and GMRES can be used in alternative. Efficiency issue is beyond the scope of our investigation but more details and comparison of performance when implementing different linear algebra solvers are found in the benchmark of the FVCA-6 Conference held in Prague, Czech Republic, in June 2011, [46].

We solve (2.1)-(2.2) on the domain for the data specified in the three test cases reported below. For all calculations, we measure the following relative errors:

error on the solution:

error on the gradient:

Test Case 1This test case is taken from [47] and is devoted to confirm that the new treatment of tensor coefficients and the standard weighted arithmetic average show the same behavior when the diffusion coefficients are regular, for example, constant or smoothly varying functions of position, even if small anisotropies along the principal direction of diffusion are present.The exact solution that we want to approximate is given by We consider two different constant diffusion tensor. The first one, called , is isotropic, while the second one, called , is anisotropic.The two diffusion tensors are given by: We solve this test case on mesh family . Each mesh of such family is built as follows. First, we remap the position of the nodes of an uniform partition by the smooth coordinate transformation: The corresponding mesh of is built from this “primal” mesh by splitting each quadrilateral cell into two triangles and connecting the barycenters of adjacent triangular cells by a straight segment. The mesh construction is completed at the boundary by connecting the barycenters of the triangular cells close to to the midpoints of the boundary edges and these latters to the boundary vertices of the “primal” mesh. The first and the second mesh of mesh family are displayed in Figure 2; mesh data are given in Table 1. Approximation errors for the isotropic diffusion tensor are shown in Figure 6; approximation errors for the anisotropic tensor are shown in Figure 7. In both plots, we also show the exact slopes proportional to and . From Figures 6 and 7, we deduce that in the case of constant diffusion tensors the performance of the new diffusion average proposed in this work and the weighted arithmetic average are almost the same. This behavior is confirmed both in the case of an isotropic diffusion tensor and in the case of an anisotropic diffusion tensor.

Table 1: Data of mesh family ; is the mesh label, is the number of cells, is the number of mesh edges, is the number of mesh nodes, and is the mesh size parameter.

1
121
400
280

2
441
1400
960

3
1681
5200
3520

4
6561
20000
13440

5
25921
78400
52480

(a)

(b)

Figure 2: First and second mesh of mesh family . This mesh family is used in test case 1.

Test Case 2In this test case, we show the behavior of the new technique that is proposed in this paper when the diffusion coefficients are discontinuous along an internal line of the computational domain. To such purpose, we split where The diffusion tensor is discontinuous across the horizontal line and is given by where and , . The exact solution is continuous across the line and is designed to ensure flux conservation, that is, continuity of the normal component of the flux field. The exact solution is given by This test case is solved on the two mesh families described below.Mesh Family This mesh sequence is the first of the mesh collection of the two-dimensional benchmark of the conference “Finite Volumes for Complex Applications—V’’ held in Aussois (France) in 2008. The first mesh and the first refined mesh of this mesh suite are shown in Figure 3; mesh data are reported in Table 2.Mesh Family Below , we consider a regular mesh formed by squares, while above we consider a regular mesh formed by squares. All the mesh nodes except those located at the boundaries and those located at the internal discontinuity line are then displaced by perturbing their position to a random position inside a square box centered at that original node position. The sides of this square box are aligned with the coordinate axis and their length is equal to (note that is the distance between two consecutive nodes in each direction). Instead, the nodes lying at are allowed to change only the abscissa in order not to modify the shape of the interface line while the nodes at the boundary are not displaced. We treat a nonmatching mesh as a conformal polygonal mesh modifying the shape of the polygons that are immediately below the line : these polygons are treated as degenerate pentagons with two parallel consecutive edges. The first mesh is built by taking and mesh refinement is implemented by doubling the value of at each refinement step, thus implying that the mesh size parameter is approximately halved when passing from one mesh to the next one in the refinement process.The first mesh and the first refined mesh of this mesh suite are shown in Figure 4; mesh data are reported in Table 3.The numerical results for the gradient and the solution approximation on mesh family are shown in the two plots of Figure 8 and on mesh family in the two plots of Figure 9. In the bottom-left corner of both plots, we show the exact slopes proportional to (gradient errors) and to and (solution errors). The solution gradient has a discontinuous normal component across such lines, and, due to this lack of regularity, the convergence rate cannot be expected to be better than that of a first-order scheme. This behavior is reflected in both left plots of Figures 8 and 9. Even if the convergence rate of the gradient is the same in the two cases, the approximation errors of the formulation using the new average are smaller than those obtained when using the standard weighted arithmetic average. Concerning the solution approximation, the results are even more spectacular because adopting the new average allows us to recover the second-order convergence rate, thus confirming the superior behavior of the new method.

Table 2: Data of mesh family ; is the mesh label, is the number of cells, is the number of mesh edges, is the number of mesh nodes, and is the mesh size parameter.

1
250
535
286

2
1000
2070
1071

3
4000
8140
4141

4
16000
32280
16281

5
64000
128560
64561

Table 3: Data of mesh family ; is the mesh label, is the number of cells, is the number of mesh edges, is the number of mesh nodes, and is the mesh size parameter.

1
56
92
37

2
224
352
129

3
896
1376
481

4
3584
5440
1857

5
14336
21632
7297

6
57344
86272
28929

7
229376
344576
115201

(a)

(b)

Figure 3: First and second mesh of mesh family . This mesh family is used in test case 2.

(a)

(b)

Figure 4: First and second mesh of mesh family . This mesh family is used in test case 2.

Test Case 3In this test case we aim at confirming the behavior of test case 2 for a more difficult case in which the diffusion tensor has intersecting discontinuities with an internal cross point.To such purpose, we split the computational domain as where The diffusion tensor is discontinuous across the horizontal line and the vertical line and is given by where and . The exact solution is continuous across such line and is designed to ensure flux conservation, that is, continuity of the normal component of the flux field. The exact solution is given by This test case is solved on the mesh family 5 that is built by displacing randomly all the nodes except those at the subdomain interfaces and of a regular partition of . We consider five meshes with ; mesh data are shown in Table 4, and the first mesh and the first refined mesh of this mesh suite are shown in Figure 5. The numerical results for the gradient and the solution approximation on mesh family are shown in the two plots of Figure 10. In the bottom-left corner of these plots, we show the exact slopes proportional to (gradient errors) and to and (solution errors). As in test case 2, the normal component of the solution gradient is discontinuous across the internal lines and . Therefore, the convergence rate that is numerically measured in this test case is proportional to , cf. the left plot of Figure 10. This fact is independent of the average technique that is used for the diffusion coefficients. Nonetheless, the approximation errors given by the new average are smaller than those obtained when using the standard weighted arithmetic average. Concerning the solution approximation, such results confirms the behavior seen in test case 2. Consequently, the finite volume approximation of the scalar solution recovers the optimal second-order convergence rate when we apply the new average technique.

Table 4: Data of mesh family ; is the mesh label, is the number of cells, is the number of mesh edges, is the number of mesh nodes, and is the mesh size parameter.

1
100
220
121

2
400
840
441

3
1600
3280
1681

4
6400
12960
6561

5
25600
51520
25921

(a)

(b)

Figure 5: First and second mesh of mesh family . This mesh family is used in test case 3.

(a)

(b)

Figure 6: Test case 1: error curves for the gradient approximation (left plot) and the solution approximation (right plot). Problem (2.1)-(2.2) is solved using mesh family and the isotropic diffusion tensor in (4.4).The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. Note that the two error curves are superimposed. In the bottom-left corners of both plots, we show the exact slopes proportional to and .

(a)

(b)

Figure 7: Test case 1: error curves for the gradient approximation (left plot) and the solution approximation (right plot). We use the anisotropic diffusion tensor of (4.4). The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. Note that the two error curves are superimposed. In the bottom-left corners of both plots, we show the exact slopes proportional to and .

(a)

(b)

Figure 8: Test case 2: error curves for the gradient approximation (left plot) and the solution approximation (right plot). Problem (2.1)-(2.2) is solved using mesh family and the discontinuous diffusion tensor (4.7).The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. In the bottom-left corner of both plots, we show the exact slopes proportional to (gradient errors) and to and (solution errors).

(a)

(b)

Figure 9: Test case 2: error curves for the gradient approximation (left plot) and the solution approximation (right plot). Problem (2.1)–(2.2) is solved using mesh family and the discontinuous diffusion tensor (4.7). The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. In the bottom-left corner of both plots we show the exact slopes proportional to (gradient errors) and to and (solution errors).

(a)

(b)

Figure 10: Test case 3: error curves for the gradient approximation (left plot) and the solution approximation (right plot). Problem (2.1)-(2.2) is solved using mesh family and the discontinuous diffusion tensor (4.11). The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. In the bottom-left corner of both plots we show the exact slopes proportional to (gradient errors) and to and (solution errors).

5. Conclusions

The numerical treatment of diffusion coefficients is an open problem and major problem, as, for example, the conductivity of different layers of soil is abruptly spatially variable. Moreover, it is particularly important in view of the problems of different scales involved and of random fields of the coefficients themselves. To this purpose, we proposed and tested a new technique for the numerical treatment of the conductivity tensor that interpolates the information coming from the different sides of a mesh interface shared by two adjacent cells. We point it out that this design is particularly suitable to the case of discontinuous conductivity tensors since the interpolation automatically adjusts its value by introducing an appropriate correction to the standard weighted average. Numerical experiments confirm this behavior. Furthermore, it may result in a particularly efficient strategy for the numerical resolution of problems where both the diffusive/dispersive and convective phenomena are simultaneously significant.

Acknowledgments

This work was funded by the PRIN 2007 project “Experimental measurement of the soil-vegetation-atmosphere interaction processes and numerical modelling of their response to climate change”. The authors also thank Dr. G. Manzini for the meshes used in test case 1 and test case 3.

References

R. Eymard, T. Gallouët, and R. Herbin, “Finite volume methods,” in Handbook of Numerical Analysis, vol. 7, pp. 713–1020, North-Holland, Amsterdam, The Netherlands, 2000.

Y. Coudière, Analyse de schémas volumes finis sur maillages non structurés pour des problèmes linéaires hyperboliques et elliptiques, Ph.D. thesis, Universitè “P. Sabatier” de Toulouse, Toulouse, France, 1999.

Y. Coudière, J.-P. Vila, and P. Villedieu, “Convergence rate of a finite volume scheme for a two-dimensional convection-diffusion problem,” Mathematical Modelling and Numerical Analysis, vol. 33, no. 3, pp. 493–516, 1999.

Y. Coudière and P. Villedieu, “Convergence rate of a finite volume scheme for the linear convection-diffusion equation on locally refined meshes,” M2AN. Mathematical Modelling and Numerical Analysis, vol. 34, no. 6, pp. 1123–1149, 2000.

E. Bertolazzi and G. Manzini, “A cell-centered second-order accurate finite volume method for convection-diffusion problems on unstructured meshes,” Mathematical Models & Methods in Applied Sciences, vol. 14, no. 8, pp. 1235–1260, 2004.

G. Manzini and A. Russo, “A finite volume method for advection-diffusion problems in convection-dominated regimes,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 13–16, pp. 1242–1261, 2008.

G. Manzini and S. Ferraris, “Mass-conservative finite volume methods on 2-D unstructured grids for the Richards' equation,” Advances in Water Resources, vol. 27, no. 12, pp. 1199–1215, 2004.

R. Haverkamp and M. Vauclin, “A note on estimating finite difference interblock hydraulic conductivity values for transient unsaturated flow problems,” Water Resources Research, vol. 15, no. 1, pp. 181–188, 1979.

N. Varado, I. Braud, P. J. Ross, and R. Haverkamp, “Assessment of an efficient numerical solution of the 1D Richards' equation on bare soil,” Journal of Hydrology, vol. 323, no. 1–4, pp. 244–257, 2006.

N. Romano, B. Brunone, and A. Santini, “Numerical analysis of one-dimensional unsaturated flow in layered soils,” Advances in Water Resources, vol. 21, no. 4, pp. 315–324, 1998.

J. M. Gastó, J. Grifoll, and Y. Cohen, “Estimation of internodal permeabilities for numerical simulation of unsaturated flows,” Water Resources Research, vol. 38, no. 12, pp. 621–6210, 2002.

A. Szymkiewicz and R. Helmig, “Comparison of conductivity averaging methods for one-dimensional unsaturated flow in layered soils,” Advances in Water Resources, vol. 34, no. 8, pp. 1012–1025, 2011.

P. Berg, “Long-term simulation of water movement in soils using mass-conserving procedures,” Advances in Water Resources, vol. 22, no. 5, pp. 419–430, 1999.

F. Chen and L. Ren, “Application of the finite difference heterogeneous multiscale method to the Richards' equation,” Water Resources Research, vol. 44, no. 7, article W07413, 2008.

L. Li, H. Zhou, and J. J. Gómez-Hernández J. Jaime, “A comparative study of three-dimensional hydraulic conductivity upscaling at the macro-dispersion experiment (MADE) site, Columbus Air Force Base, Mississippi (USA),” Journal of Hydrology, vol. 404, no. 3-4, pp. 278–293, 2011.

W. E and B. Engquist, “The heterogeneous multiscale methods,” Communications in Mathematical Sciences, vol. 1, no. 1, pp. 87–132, 2003.

Y. Efendiev and T. Y. Hou, Multiscale Finite Element Methods. Theory and Applications, vol. 4, Springer, New York, NY, USA, 2009.

I. Babuška, “Solution of interface problems by homogenization—I, II,” SIAM Journal on Mathematical Analysis, vol. 7, no. 5, pp. 603–634, 1976.

I. Babuška, “Homogenization and its application. Mathematical and computational problems,” in Numerical Solution of Partial Differential Equations III, B. Hubbard, Ed., pp. 89–116, Academic Press, New York, NY, USA, 1976.

I. Babuška, “Solution of interface problems by homogenization—III,” SIAM Journal on Mathematical Analysis, vol. 8, no. 6, pp. 923–937, 1977.

L. J. Durlofsky, Y. Efendiev, and V. Ginting, “An adaptive localglobal multiscale finite volume element method for two-phase flow simulations,” Advances in Water Resources, vol. 30, pp. 576–588, 2007.

P. Dostert, Y. Efendiev, and B. Mohanty, “Efficient uncertainty quantification techniques for Richards’ equation,” Advances in Water Resources, vol. 32, pp. 329–339, 2009.

L. Bergamaschi, S. Mantica, and G. Manzini, “A mixed finite element—finite volume formulation of the black-oil model,” SIAM Journal on Scientific Computing, vol. 20, no. 3, pp. 970–997, 1999.

C. Gallo and G. Manzini, “2-D numerical modeling of bioremediation in heterogeneous saturated soils,” Transport in Porous Media, vol. 31, no. 1, pp. 67–88, 1998.

C. Gallo and G. Manzini, “A mixed finite element finite volume approach for solving biodegradation transport in groundwater,” International Journal for Numerical Methods in Fluids, vol. 26, no. 5, pp. 533–556, 1998.

C. Gallo and G. Manzini, “A fully coupled numerical model for two-phase flow with contaminant transport biodegradation kinetics,” Communications in Numerical Methods in Engineering, vol. 17, no. 5, pp. 325–336, 2001.

L. Beirão da Veiga, V. Gyrya, K. Lipnikov, and G. Manzini, “Mimetic finite difference method for the Stokes problem on polygonal meshes,” Journal of Computational Physics, vol. 228, no. 19, pp. 7215–7232, 2009.

A. Bouziani and I. Mounir, “On the solvability of a class of reaction-diffusion systems,” Mathematical Problems in Engineering, vol. 2006, Article ID 47061, 15 pages, 2006.

A. Cangiani and G. Manzini, “Flux reconstruction and solution post-processing in mimetic finite difference methods,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 9–12, pp. 933–945, 2008.

W. Correa de Oliveira Casaca and M. Boaventura, “A decomposition and noise removal method combining diffusion equation and wave atoms for textured images,” Mathematical Problems in Engineering, vol. 2010, Article ID 764639, 20 pages, 2010.

M. Dehghan, “On the numerical solution of the diffusion equation with a nonlocal boundary condition,” Mathematical Problems in Engineering, vol. 2003, no. 2, pp. 81–92, 2003.

S. Hu, Z. Liao, D. Sun, and W. Chen, “A numerical method for preserving curve edges in nonlinear anisotropic smoothing,” Mathematical Problems in Engineering, vol. 2011, Article ID 186507, 14 pages, 2011.

Z. Liao, S. Hu, D. Sun, and W. Chen, “Enclosed Laplacian operator of nonlinear anisotropic diffusion to preserve singularities and delete isolated points in image smoothing,” Mathematical Problems in Engineering, vol. 2011, Article ID 749456, 15 pages, 2011.

R. J. Moitsheki, “Transient heat diffusion with temperature-dependent conductivity and time-dependent heat transfer coefficient,” Mathematical Problems in Engineering, vol. 2008, Article ID 347568, 9 pages, 2008.

N. P. Moshkin and D. Yambangwai, “On numerical solution of the incompressible Navier-Stokes equations with static or total pressure specified on boundaries,” Mathematical Problems in Engineering, vol. 2009, Article ID 372703, 26 pages, 2009.

W. Wenquan, Z. Lixiang, Y. Yan, and G. Yakun, “Finite element analysis of turbulent flows using les and dynamic subgridscale models in complex geometries,” Mathematical Problems in Engineering, vol. 2011, Article ID 712372, 20 pages, 2011.

J. Yu, Y. Yang, and A. Campo, “Approximate solution of the nonlinear heat conduction equation in a semi-infinite domain,” Mathematical Problems in Engineering, vol. 2010, Article ID 421657, 24 pages, 2010.

Y. Zheng and Z. Zhao, “A fully discrete Galerkin method for a nonlinear space-fractional diffusion equation,” Mathematical Problems in Engineering, vol. 2011, Article ID 171620, 20 pages, 2011.

P. Grisvard, Elliptic Problems in Nonsmooth Domains, vol. 24, Pitman, Boston, Mass, USA, 1985.

E. Bertolazzi and G. Manzini, “A second-order maximum principle preserving finite volume method for steady convection-diffusion problems,” SIAM Journal on Numerical Analysis, vol. 43, no. 5, pp. 2172–2199, 2005.

E. Bertolazzi and G. Manzini, “On vertex reconstructions for cell-centered finite volume approximations of 2D anisotropic diffusion problems,” Mathematical Models & Methods in Applied Sciences, vol. 17, no. 1, pp. 1–32, 2007.

G. Manzini and M. Putti, “Mesh locking effects in the finite volume solution of 2-D anisotropic diffusion equations,” Journal of Computational Physics, vol. 220, no. 2, pp. 751–771, 2007.

G. P. Leonardi, F. Paronetto, and M. Putti, “Effective anisotropy tensor for the numerical solution of flow problems in heterogeneous porous media,” in Proceedings of the CMWR 16th International Conference, Copenhagen, Denmark, 2006.

M. Quintard and S. Whitaker, “Single-phase flow in porous media: effect of local heterogeneities,” Journal de Mecanique Theorique et Appliquee, vol. 6, no. 5, pp. 691–726, 1987.

L. Tartar, “Estimations fines des coefficients homogénéisés,” in Ennio De Giorgi Colloquium, vol. 125, pp. 168–187, Pitman, Boston, Mass, USA, 1985.

R. Eymard, G. Henri, R. Herbin, R. Klofkorn, and G. Manzini, “3D Benchmark on discretizations schemes for anisotropic diffusion problems on general grids,” in Finite Volumes for Complex Applications VI, Problems and Perspectives, J. Fort, J. Furst, J. Halama, R. Herbin, and F. Hubert, Eds., vol. 2, pp. 95–130, Springer, 2011.

L. Beirão da Veiga and G. Manzini, “A higher-order formulation of the mimetic finite difference method,” SIAM Journal on Scientific Computing, vol. 31, no. 1, pp. 732–760, 2008.