Comparison of DDFV and DG Methods for Flow in Anisotropic Heterogeneous Porous Media

We consider unsteady single-phase flows modeled by the Richards' equation. We compare Discrete Duality Finite Volume and Discontinuous Galerkin schemes applied to discretize the diffusive term. Backward Differentiation Formula of second order is used for the time stepping method. Accuracy and robustness of these two types of schemes are tested and compared on various test cases, especially in anisotropic heterogeneous media.


INTRODUCTION
Prediction of fluid flows in geological subsurface is a key process in many applications like groundwater hydrology, oil recovery, and more recently CO 2 geological storage and deep geothermal science. Several works are devoted to design numerical schemes (listed below) for simulating these types of flows. This article is a comparison of two conservative methods which are accurate (i.e: having a second-order convergence), robust in anisotropic heterogeneous media (i.e: with a diffusivity tensor inducing a space-dependent preferential direction for the flow) and valid for general meshes (i.e: not only for K-orthogonal meshes). We present a preliminary work dedicated to single-phase flows. An extension of this study is to consider two phases (gas and water) in the frame of CO 2 injection in deep aquifers.
A various choice of methods is available to approximate nonlinear diffusion equations. For instance to discretize the Richards equation, Celia et al. (1990) consider Finite Element (FE), Narasimhan and Witherspoon (1976) apply Finite Volume (FV) with two-point flux approximation, Manzini and Ferraris (2004) use diamond cell FV, and Knabner and Schneid (2002) employ Mixed Finite Element (MFE). Several other FV methods and variants exist, including the MPFA, SUSHI and VAG schemes, see Eymard et al. (2012) and references therein for details. In this work, we propose an alternative approach by applying and comparing DDFV and DG methods. The DDFV methodology was introduced about ten years ago by Hermeline (2000). It has been applied to a variety of problems over the past few years, from isotropic scalar diffusion with Domelevo and Omnes (2005) to electrocardiology, Coudie`re et al.
(2009), among others. This FV method is a natural choice because it has proven robustness when applied to anisotropic and heterogeneous diffusion problems on general, most notably non-matching and distorted meshes. It is locally conservative and preserves the symmetry of the continuous problem, which allows efficient iterative solvers to be employed. Besides it offers a L 2 -norm superconvergence on various examples, together with a very accurate appoximation of the gradient, as was noted in Herbin and Hubert (2008). In this paper, we adapt a version of the DDFV scheme described in Coudie`re et al. (2009), Krell (2010) which allows discontinuities of the permeability tensor along faces. DG methods have appeared thirty years ago with Lesaint and Raviart (1974) and are efficient to approximate equations present in engineering sciences. Advantages of these methods are multiple especially the local conservation (like FV and MFE). The order of accuracy (like FE and MFE) can be different for each element of the mesh, facilitating the p-refinement. The flexibility in the use of non-matching meshes (like FV) is particularly adaptated for the h-refinement. Several DG methods can be used for one-phase and two-phase flows in porous media, such as the local DG method described in Bastian et al. (2007), Fagherazzi et al. (2004) and the nonsymmetric or the symmetric interior penalty Galerkin method used in Bastian (2003), Klieber andRivie`re (2006), Sochala et al. (2009). Here, we choose the Symmetric Weighted Interior Penalty DG method (SWIP), because it preserves the natural symmetry in the discrete diffusion operator and uses diffusion-dependent weighted averages to reduce the amount of stabilization required, Di Pietro and Ern (2012). Regarding time discretization, we propose to use the second-order backward differentiation formula (BDF2) from Curtiss and Hirschfelder (1952). This implicit scheme avoids the CFL condition particularly restrictive for explicit schemes when the diffusion term is present. Furthermore, BDF2 is more efficient than implicit Runge-Kutta (DIRK) schemes, which require several stages (three for the second-order version). Besides, the last stage of DIRK schemes is ill-posed where the soil is saturated, as was underlined in Sochala et al. (2009). Finally, BDF2 is preferred to Crank-Nicolson which is not L-stable, Hairer and Wanner (2010).
The aim of this work is to test and compare the DDFV and SWIP schemes for solving flows in anisotropic heterogeneous porous media. We propose to consider the conservative form of the Richards equation in a two-dimensional domain X: where T is the final time of simulation, @X D (resp. @X N Þ the boundary of X where a Dirichlet condition w D (resp. Neumann condition v N ) is applied, and n X the outward unit normal of @X ¼ @X D [ @X N . The vector e z ¼ ð0; 1Þ t is the gravity contribution. The pressure head w is the unknown and is related to the velocity v through Darcy's law vðwÞ ¼ ÀKðwÞðrw þ e z Þ. Two constitutive relationships w 7 ! hðwÞ and w 7 ! KðwÞ are necessary to close the model. Examples of water content h and conductivity K are indicated in Sections 3.1 and 3.2. For the sake of simplicity, we suppose that the source term f is equal to zero (except for the validation test case described in Sect. 3.1). We investigate three Test Cases (TC) to compare DDFV and SWIP methods in several configurations. TC1 and TC2 are column infiltration problems in an isotropic homogeneous soil. TC1 has an analytical solution and permits to validate the convergence rate of the methods. The stencil of the DDFV method produces a sparser matrix than SWIP. In return, the latter leads to a narrower bandwith and smaller condition numbers. TC2 features a stiff pressure front which highlights significant differences for coarse meshes. DDFV underestimates the front propagation speed while SWIP presents nonphysical oscillations. TC3 is a quarter five-spot configuration in an anisotropic heterogenous soil. The two methods only differ in the evaluation of the mass error, depending on the choice of the nonlinear solver tolerance. The outline is as follows. Section 1 introduces the space discretization associated with DDFV and SWIP methods. Section 2 describes the time discretization performed by the BDF2 formula and the Picard fixed point algorithm used to solve the nonlinear discrete system. Section 3 details the results on the three TC described above. Conclusions and perspectives are drawn after.

Notations
Let fT h g h>0 be a shape-regular family of unstructured meshes of X consisting for simplicity of affine triangles. For an element s 2 T h , let @s denote its boundary and n s its outward unit normal. The set F h of mesh faces is is the set of Dirichlet (resp. Neumann) faces. For a face r 2 F i h , there exist s þ and s À in T h such that r ¼ @s þ \ @s À . We denote by F s & F h the set of the faces r such that @s ¼ [ r2F s r, and by F s & F h the set of the faces r that share s as a vertex. We define n r as the unit normal to r pointing from s À to s þ (Fig. 1).
The set S h of mesh vertices is partitioned into gathers all the endpoints of the edges in F N h and S D h are the remaining boundary vertices. A secondary mesh is associated to this set S h as follows: any vertex s 2 S i h [ S N h is associated to a unique polygonal cell S whose vertices are the centers of the triangles s and the midpoints of the faces r that share s as a vertex; and which contains s (Fig. 2). Consequently, each face r is associated to its two neighbouring cells s þ and s À and its two vertices, denoted by s þ and s À . The following orientation assumption is made for these notations: although the complete orientation is arbitrary and irrelevant in the sequel. Here, x s is the center of gravity of s 2 T h and x s are the coordinates of s 2 S h . We also denote by x r the midpoints of the faces r 2 F h (Fig. 3). Note that a face r & @X has only one neighbouring triangle which is denoted by s À in the sequel, so that the normal n r is also the unit normal to @X outward of X.
The DDFV method couples, through the flux computations, two finite volume schemes written on a primary mesh T h (considering the cells s as control volumes) and on a secondary mesh S h (considering the cells S as control volumes). Hence, the DDFV discrete unknown is a set of degrees of freedom W h , that eventually defines a function w h (Sect. 1.2). The DG method directly defines a discrete unknown function w h (Sect. 1.3) and is based on a discrete variational formulation.

Discrete Duality Finite Volume Method
The DDFV method approximates the average of w on each control volume s 2 T h and on each control volume S associated to a vertex s 2 S i h [ S N h . Hence, we define the discrete DDFV unknown as: Several functions can be associated to these degrees of freedom, for instance the usual piecewise constant function on each s 2 T h (or its sibling piecewise constant on each S for s 2 S h ), but also the function w h that is piecewise linear on the triangles with base r and vertex x s (for all s 2 T h and r 2 F s ) and such that w h ðx s Þ ¼ w s and w h ðx s Þ ¼ w s .
If we integrate Equation (1) on a control volume s 2 T h we get: Note that the flux ÀKðwÞðrw þ e z Þ Á n s is continuous through each face r of s. The DDFV scheme reads, for any control volume s 2 T h : where w r;s (resp. r r;s W h ) approximates w (resp. rw) on the side s of the face r, the normal N r ¼ jrjn r is the normal to r outward of s accounting for the length of r and V r approximates R r v N . Note that in the TC presented in Section 3, the boundary data is piecewise constant and we can take exactly V r ¼ R r v N . There are several possible choices for the value w r;s . We take here: because it seems to provide the best accuracy and robustness. For each vertex s 2 S D h , we choose to take w s ¼ w D ðx s Þ from the Dirichlet boundary data. For any control volume S associated to a vertex s 2 In these summations, if r 2 F Ã s then the vertex s is identified to s À , the neighbouring triangles are s AE and the interface between the control volumes associated to s À and s þ has two parts e þ and e À , of which the corresponding normals are N e AE ¼ n e AE je AE j as depicted in Figure 3. If r 2 F N s , then there is only one triangle s neighbouring r but the flux still has two terms, one along the edge e defined as in the previous case (with the normal N e ¼ n e jej) and one along the line ðx s ; x r Þ & r. The number V rs approximates the integral of the Neumann boundary data along this line. Given a face r 2 F i h , we denote by: the fluxes out of s À and s þ through the face r and by: is the flux from s to the opposite endpoint of r. The flux F e is continuous through the interface between S À and S þ by construction. We ensure the continuity of the fluxes F r;s AE through r by introducing an auxiliary unknown w r used to compute the gradients r r;s À W h and r r;s þ W h , following the DDFV strategy from Coudie`re et al. (2009): Notations for DDFV method.
where 2D r;s AE ¼ AE detðx s AE À x r ; x s þ À x s À Þ, meaning that D r;s AE are the areas of the triangles with base r and vertices x s AE . In the boundary case r 2 F N h , there is only one gradient, namely r r;s À W h . The auxiliary unknown is then the unique solution to the linear equation F r;s À þ F r;s þ ¼ 0 in the general interior case and to the linear equation F r;s À ¼ V r in the boundary case. Then the flux between the primary cells s AE is F r ¼ F r;s À ¼ ÀF r;s þ and the scheme (2, 3) reads: In view of the expression of the fluxes and the gradients, we adopt some additional notations in the general interior case. On each side s AE of r, we define the following Gram matrices: and a AE , b AE are defined similarly by: and can be easily solved. Afterwards, we replace w r by its expression and find that the calculation of the fluxes F r and F e can be carried out as follows: In the case of a face r 2 F N h \ s, we assume that the neighbouring triangle s is s À and the equation on w r reads a À ðw r À w s À Þ þ b À ðw s þ À w s À Þ þ a À ¼ V r . It is solved easily and results in the expression: Lastly, if r 2 F D h , the quantity F e does not need to be defined, as only F r adds a contribution to the system. In this case, the auxiliary w r is chosen as w r ¼ 1 2 ðw s þ þ w s À Þ. By assembling these local contributions, the semidiscrete DDFV scheme finally reads: , M h is a diagonal mass matrix and A h ðW h Þ is a symmetric and positive matrix.

Discontinuous Galerkin Method
The Discontinuous Galerkin methods are based on variational formulations (as for the classical finite element). The discontinuous finite element space V h is defined as: where P p ðsÞ is the set of polynomials of degree less than or equal to p on an element s. We observe that the functions in V h are not necessarily continuous. This fact is exploited by selecting basis functions which are locally supported in a single mesh element. As a consequence, when a triangular mesh and a P 1 -Lagrange approximation are chosen, the unknowns are localized at the vertices of each triangle as illustrated in Figure 4. A linear approximation yields the same order of convergence as DDFV, though higher-order elements could also be considered. We define the weighted average operator {} x,r and the jump operator [ ] r as follows: for a function which is possibly two-valued on r: fng x;r ¼ def x À n À þ x þ n þ and ½n r ¼ def n À À n þ where n AE ¼ nj s AE . For boundary faces, fng x;r ¼ nj r and ½n r ¼ nj r . For vector-valued functions, average and jump operators are defined componentwise. For the Richards equation, the SWIP method can be written as: A specific choice of the weights is made to yield robust error estimates with respect to the diffusivity: In the penalty term, g is a positive parameter (to be taken larger than a minimal threshold depending on the shape-regularity of T h ), and d r is the diameter of the face r, i.e. the largest diameter of the triangle(s) of which r is a face. The coefficient c K is defined as: with d Kn ¼ Kn r Á n r . The above choice of x À and x þ leads to the harmonic mean of the diffusivity tensor: Using harmonic means is particularly adapted at an interface between poorly and highly conductive media. In this case, the harmonic mean value of the diffusivity is close to the value in the poorly conductive medium, Di Pietro and Ern (2012).
Again, by assembling the local contributions, the semi-discrete DG scheme takes the form of (4), where M h is a block diagonal matrix, A h is a block matrix and W h is the set of the values of w h at each degree of freedom.

Backward Differentiation Formula
The next step consists in approximating the time derivative in Equation (4). We propose the second order Backward Differentiation Formula (BDF2): where the superscript n denotes the time level, and dt is the constant time step, chosen such that T =dt is an integer. Substituting the last expression in Equation (4) finally yields the following fully discrete scheme: The vector W 1 h is computed by the Crank-Nicolson scheme.

Linearization Algorithm
The above is obviously a nonlinear system, and may thus be solved through an iterative procedure, which we Unknowns localization for each method.
describe in this section. The iteration level will be reffered to as the superscript m. As suggested in Manzini and Ferraris (2004) We use the LU decomposition for solving each linear system and our stopping criterion is jjdW n;m h jj 2 jjW nÀ1 h jj À1 2 , where is the user-defined tolerance.

Mass Conservation Properties
The study of the mass evolution and repartition (i.e. the separation over time of the mass variation in the domain, the mass inflow and the mass outflow) helps to understand the behavior of the system. It also allows to quantify the mass defect generated by the linearization algorithm. The volume of water in the domain X at time ndt is obtained by integrating the volumetric water content in X: By using the BDF2 formula, we have: where the quantities F n D and F n N are defined as: and n stands for the numerical error in the resolution of the nonlinear system. The reconstructed normal velocity v Ã;n h on a Dirichlet face r is estimated from the pressure: v Ã;n h j r ¼ ÀKðw n r;s À Þ r r;s À W n h þ e z À Á Á n r for DDFV vðw n h j r Þ Á n r þ gc K d r ðw n h j r À w D Þ for SWIP ( The variation of the volume of water over the time step ½ðn À 1Þdt; ndt is deduced from a reformulation of Equation (5), yielding: where the flux U n 2 fU n D ; U n N g and the error e n over the time step ½ðn À 1Þdt; ndt are defined as: The initialization of the fluxes recursive formulas depends on the scheme used at the first time step. For the Crank-Nicolson scheme, the global mass conservation is written as:

RESULTS
We propose three TC with increasing difficulties corresponding to a more general form of the conductivity defined as: where kðwÞ is the relative permeability, KðxÞ the intrinsic tensor permeability of the soil; q, g and m are respectively the fluid density, the gravity constant and the dynamic viscosity. We first consider two downward infiltrations in an isotropic homogeneous medium (i.e. K ¼ I). The first TC presents an analytical solution to verify the theoretical convergence rates and the second TC is the propagation of a stiff pressure front. We also study the quarter five-spot problem in an anisotropic heterogeneous medium (i.e. K 6 ¼ I and is space-dependent).

Isotropic Homogeneous Validation Test Case
We propose a TC with an analytical solution to compare the two methods in terms of matrix properties and convergence. The domain is X ¼ ½0; 4 Â ½0; 20 (in cm) and the final time T ¼ 2 min. The analytical solution: wðz; tÞ ¼ 20:4 tanh½0:5ðz þ t 12 À 15Þ À 41:1 is used to determine adequate source term and Dirichlet boundary condition enforced on @X. The water content and the conductivity are defined by Haverkamp's constitutive relationships which are derived in Haverkamp et al. (1977): K s ¼ 9:44 Â 10 À3 cm:s À1 ;Ã ¼ 0:0524 cm À1 ; c ¼ 4:74: For a triangular mesh, the number of unknowns N u is: where N t is the number of triangles, N n the number of nodes and N D n the number of nodes where a Dirichlet condition is enforced (and which is equal to the number of nodes located on @X). Consequently, the number of unknowns for DDFV is about half the number of unknowns for SWIP when the mesh is sufficiently fine (i.e. N n ) N D n ): because N t ' 2N n from the Euler relations. We have therefore constructed two families of unstructured meshes fM i g 16i66 and fM i g 16i66 respectively for DDFV and SWIP, and such that the number of unknowns for M i and M i are close. For each mesh, Table 1 refers to N n ; N t and N u . The number of nonzero elements nnz, the mean number of nonzero elements per row nnz (¼ nnz=N u Þ and the bandwidth bw of the corresponding matrix (obtained after a reverse Cuthill-McKee ordering, Cuthill and McKee (1969)) are also indicated. The first remark is that the number of nonzero elements is different for each method since it depends on the stencil. By considering fine meshes (typically M 5 , M 6 , M 5 and M 6 ), we have for DDFV: while we have for SWIP: The second remark is that the bandwidth is smaller for SWIP which has a more simple connectivity graph due to the block structure of its matrix. To sum up, DDFV requires less storage capacity whereas SWIP seems more efficient for solving each linear system. Table 2 presents the errors on the hydraulic head e w;X and on the velocity e v;X for all the meshes: ; e v;X ¼ jjv À v h jj L 2 ðXÂ½0;T Þ jjvjj L 2 ðXÂ½0;TÞ For the two methods, the results confirm theoretical estimations since a second-order convergence is observed on the hydraulic head and a first-order convergence is verified on the velocity. Table 3 reports the mean number N it of iterations in the linearization algorithm and the mean condition number j 1 2 : where N n it is the number of iterations performed at time level n and k n;m max (resp. k n;m min ) is the maximal (resp. minimal) eigenvalue of the matrix at the m-th nonlinear iteration of the n-th time step. We observe that the mean condition number is inversely proportional to the mesh size h. This is consistent with a finite difference result which states that the condition number j 2 is proportional to h À2 . Furthermore, the SWIP method significantly reduces e w;X (by a factor of 1.5) and j 1 2 (by a factor of 6) obtained with the most classical symmetric interior penalty Galerkin method (which uses x þ ¼ x À ¼ 0:5). To balance the various error sources (due to space and time discretizations as well as nonlinear systems resolution), an adaptive inexact Newton method and adaptive time-stepping based on a posteriori estimation could be implemented, as was done in Ern and Vohralı´k (2010).

Isotropic Homogeneous Stiff Case
This TC is a stiff infiltration problem proposed by Celia et al. (1990). The domain is X ¼ ½0; 20 Â ½0; 100 (in cm) and the final time is T ¼ 48 h. A constant initial condition is considered. A Dirichlet condition is imposed on the top T and the bottom B of the column. An homogeneous Neumann condition (corresponding to a zero flux) is imposed on the lateral parts L of the domain (Fig. 5):  The water content and the conductivity are defined by Van Genuchten's constitutive relationships (van Genuchten, 1980), and plotted in Figure 6: with parameters: h s ¼ 0:368; h r ¼ 0:102; K s ¼ 9:22 Â 10 À3 cm:s À1 ; The evoked ''stiffness'' of this TC is related to the strong constrained overpressure (equal to 9:25 m) imposed on the top of the column. A sharp variation of the conductivity ðKðÀ75 cmÞ=KðÀ10 mÞ ¼ 8:92 Â 10 4 Þ is induced, Figure 6.
We focus here on the horizontal mean value of the pressure head w h z ð Þ: 8z 2 ½0; 100; w h ðzÞ ¼ 1 20 . The DDFV pressure profiles are depicted by solid lines and the SWIP pressure profiles are depicted by dashed lines. Each method presents one drawback when coarse meshes are used. For DDFV, a delay of the pressure front is observed, especially when M 4 is employed since a significant shift is observed between the two methods. For SWIP, non-physical oscillations appear and could be reduced with the help of a slope limiter as in Cockburn and Shu (1998). In this stiff case, only DDFV satisfies the discrete maximum principle whereas SWIP has a better estimation of the propagation speed of the pressure front. When the meshes are sufficiently fine, the two methods yield the same pressure in agreement with Celia et al. (1990) and Manzini and Ferraris (2004).

Anisotropic Heterogeneous Five-Spot Problem
This TC is inspired by the quarter five-spot problem studied in Simmons et al. (1959), which reproduces an elementary cell of a periodic network consisting of sources and sinks. The domain X ¼ ½0; 1 2 (in m) is divided into four parts (Fig. 8): The soil properties are the same as those of the previous TC except for the piecewise constant intrinsic permeability: where 1 X i ðxÞ is the indicator function of the set X i , D is a diagonal matrix and R x i is the rotation matrix associated to X i : The angles of rotation are: The matrix R t x denotes the transpose of R x . This permeability induces a preferential direction for the flow on each part of the domain as illustrated in Figure 8.
The final time is T ¼ 6 h. An hydrostatic initial condition is considered, an homogeneous Neumann condition is imposed on the boundary of the domain except on the top right corner C where an incoming flux v N is enforced TC2 -pressure head w h ðzÞ at 24 h and 48 h obtained with different meshes for DDFV (solid line) and SWIP (dashed line). and on the line L ¼ ½0; 0:025 Â f0g where a zero pressure condition is applied: The condition v N ¼ À5 Â 10 À3 cm:s À1 on C corresponds to a water injection of 9 kg per hour. Figure 9 shows that the triangles never cross the interface between each subdomain. The meshes are obtained from the FreeFrem Ó software, which allows to define anisotropic heterogeneous metrics. Figure 10 plots isolines of the overpressure w T h À w 0 h at the final time of the simulation. The two methods yield the same results for fine meshes (8 760 triangles for DG, 16 446 for DDFV). Concerning the convergence of the solution as h tends to 0, SWIP needs further refinement than DDFV when isotropic meshes are used. Therefore, DDFV is less sensitive than SWIP to the choice of the mesh when an anisotropic permeability is considered. Figure 11 and Figure 12 present results on mass conservation. Multiplying Equations (6) and (7) by the water density q, and summing over the time intervals in ½0; ndt leads to: ÁM n is the total mass variation over the time interval ½0; ndt and the quantities P n i¼1 M i in , P n i¼1 M i out and P n i¼1 jE i j are respectively the total water inflow, the total water outflow and the total mass balance defect cumulated at time ndt. The quantities ÁM n ; P n i¼1 M i in and P n i¼1 M n out are close for the two methods and three phases are clearly recognizable in Figure 11. The first phase ½0; 1 h corresponds to the soil saturation: the total mass variation equals the total water inflow. The second phase ½1 h; 4 h features soil saturation and drainage : the total mass variation increases more slowly than during the previous phase because exfiltration occurs. The last phase ½4 h; 6 h is the soil drainage only since the injection is stopped, inducing a diminution of the total mass variation.
The total mass balance defect (or mass error) P n i¼1 jE i j is plotted in Figure 12, where various tolerances in the linearization algorithm are tested. When ¼ 10 À3 , the defects look similar, but lowering the tolerance has no improving effect with SWIP. Mass error can be reduced by using higher order approximations instead. Meanwhile, about two orders of magnitude can still be gained with DDFV before convergence, which is achieved for ¼ 10 À5 .

CONCLUSION
We presented a comparison between a discontinuous finite element method and a more recent finite volume scheme on various test cases. As was presented in Section 1, the approaches are quite different and the numerical results highlight different behaviors. First, DDFV provides a sparser structure than SWIP which produces in return a narrower bandwith and a better condition number. In the stiff case, DDFV verifies the discrete maximum principle while SWIP evaluates better the speed of propagation for coarse meshes. Finally, in the quarter five-spot problem, the two methods are in accordance for the overpressure and the mass repartition over time. The mass balance defect can be considerably lowered with DDFV by requiring a more restrictive tolerance in the Picard algorithm.