We develop a mixed finite element method for single phase flow in porous media that reduces to cell-centered finite differences on quadrilateral and simplicial grids and performs well for discontinuous full tensor coefficients. Motivated by the multipoint flux approximation method where sub-edge fluxes are introduced, we consider the lowest order Brezzi-Douglas-Marini (BDM) mixed finite element method. A special quadrature rule is employed that allows for local velocity elimination and leads to

more »

... a symmetric and positive definite cell-centered system for the pressures. Theoretical and numerical results indicate second-order convergence for pressures at the cell centers and first-order convergence for sub-edge fluxes. Second-order convergence for edge fluxes is also observed computationally if the grids are sufficiently regular. Keywords: mixed finite element, multipoint flux approximation, cell centered finite difference, tensor coefficient, error estimates AMS Subject Classification: 65N06, 65N12, 65N15, 65N30, 76S05 1. Introduction. Mixed finite element (MFE) methods have been widely used for modeling flow in porous media due to their local mass conservation and accurate approximation of the velocity. They also handle well discontinuous coefficients. A computational drawback of these methods is the need to solve an algebraic system of saddle point type. One possible approach to address this issue is to use the hybrid form of the MFE method [9, 15] . In this case the method can be reduced to a symmetric positive definite system for the pressure Lagrange multipliers on the element faces. Alternatively, it was established in [29] that, in the case of diagonal tensor coefficients and rectangular grids, MFE methods can be reduced to cell-centered finite differences (CCFD) for the pressure through the use of a quadrature rule for the velocity mass matrix. This relationship was explored in [33] to obtain convergence of CCFD on rectangular grids. This result was extended to full tensor coefficients and logically rectangular grids in [7, 6] , where the expanded mixed finite element (EMFE) method was introduced. The EMFE method is very accurate for smooth grids and coefficients, but loses accuracy near discontinuities. This is due to the arithmetic averaging of discontinuous coefficients. Higher order accuracy can be recovered if pressure Lagrange multipliers are introduced along discontinuous interfaces [6] , but then the cell-centered structure is lost. Several other methods have been introduced that handle well rough grids and coefficients. The control volume mixed finite element (CVMFE) method [16] is based on discretizing the Darcy's law on specially constructed control volumes. Mimetic finite difference (MFD) methods [23] are designed to mimic on the discrete level critical properties of the differential operators. The approximating spaces in both methods are closely related to RT 0 , the lowest order Raviart-Thomas MFE spaces [27] . These relationships have been explored in [17, 30] and [10, 12] to establish convergence of the CVMFE methods and the MFD methods, respectively. However, as in the case of MFE methods, both methods lead to an algebraic saddle point problem. The multipoint flux approximation (MPFA) method [1, 2, 19, 20] has been developed as a finite volume method and combines the advantages of the above mentioned methods, i.e., it is accurate for rough grids and coefficients and reduces to a cell-centered stencil for the pressures. However, due to the non-variational formulation of the MPFA, there