Developing a predictive understanding of the terrestrial water cycle at local to global scale is essential for accurate assessment of water resources. Higher spatial resolution in the land component of the Energy Exascale Earth System Model alone is insufficient to meet the U.S. Department of Energy's 10-year vision. Next-generation terrestrial models need to move beyond one-dimensional systems by including scale appropriate three-dimensional physics formulations. To this end, we are developing a terrestrial dynamical core which is tailored to solving thermal mass balance equation at global scales. The method of choice should converge quadratically in the pressure and velocity in the presence of discontinuous coefficients and non-orthogonal grids.

Given these requirements, we study strong scaling of the solution of linear systems resulting from various choices of discretization method of the Poisson equation in mixed form. While mixed finite element methods are a natural choice, they lead to a saddle point system which presents computational challenges. We compare the lowest order Brezzi-Douglas-Marini space (BDM) against a method by Wheeler and Yotov (2006) [M. F. Wheeler and I. Yotov, A multipoint flux mixed finite element method, SIAM J. Numer. Anal. 44:5 (2006) 2082-2106.]. In this method, the drawback of the resulting saddle point problem is overcome by using a vertex quadrature rule, decoupling the velocity systems around the vertices of the mesh. This allows for the velocity to be eliminated locally, resulting in a symmetric positive definite system in terms of pressure only.