Simulation of flow in porous media requires the numerical approximation of elliptic partial differential equations. Mixed finite element methods are frequently employed, because of local mass conservation and accurate approximation of both pressure and velocity. Mixed methods give rise to "elliptic" saddle-point (ESP) matrices, which are difficult to solve numerically. In addition, the problems to be modelled in ground water flow require that the hydraulic conductivity or absolute permeability be a tensor, which adds additional complexity to the resulting saddle-point matrices. This research develops several preconditioners for restarted GMRES solution of the ESP linear systems. These preconditioners are based on a new class of polynomial matrices, which we refer to as Multi-band Jacobi Polynomial (JMP) matrices. Applications of these preconditioners to the numerical solution of two and three spatial dimensional flow equations with tensor coefficients using rectangular lowest order Raviart-Thomas spaces are presented.