We consider a covolume method for a system of first order PDEs resulting from the mixed formulation of the variable-coefficient-matrix Poisson equation with the Neumann boundary condition. The system may be used to represent the Darcy law and the mass conservation law in anisotropic porous media flow. The velocity and pressure are approximated by the lowest order Raviart-Thomas space on rectangles. The method was introduced by Russell [Rigorous Block- centered Discretizations on Irregular Grids: Improved Simulation of Complex Reservoir Systems, Reservoir Simulation Research Corporation, Denver, CO, 1995] as a control-volume mixed method and has been extensively tested by Jones [A Mixed Finite Volume Elementary Method for Accurate Computation of Fluid Velocities in Porous Media, University of Colorado at Denver, 1995] and Cai et al. [Computational Geosciences, 1 (1997), pp. 289-345]. We reformulate it as a covolume method and prove its first order optimal rate of convergence for the approximate velocities as well as for the approximate pressures
