Mixed Finite Element Methods (FEM) and multipoint flux Finite Volume Methods (FVM) can be used to accurately discretize Richards equation for porous media with heterogeneous coefficients and distorted meshes. Techniques such as inexact quadrature (lumping) can be used to make more solver-friendly and efficient formulations in exchange for some problem-dependent accuracy. Different multilevel solvers such as balancing domain decomposition by constraints and algebraic multigrid are suitable for these different formulations. This work investigates tradeoffs in accuracy versus cost metrics such as number of degrees of freedom, memory usage/bandwidth requirements, and solver efficiency using the best available scalable techniques. We consider a range of test problems with varying degrees of heterogeneity and compare accuracy metrics relevant to hydrologic applications. Results using the multipoint flux FVM and mixed FEM with exact and inexact quadrature are analyzed using the method of manufactured solutions and heterogeneous benchmark problems.