Groundwater, which accounts for 30% of freshwater reserves globally, is a vital source for human water supply. Climate change is expected to impact the quality and quantity of groundwater in the future. Numerous observational and modeling studies have shown a positive soil moisture-rainfall feedback. Spatial variability of topography, soils, and vegetation play a significant role in determining the response of land surface states (soil moisture) and fluxes (runoff, evapotranspirtiaon). Despite the obvious need to accurately represent soil moisture dynamics, the current version of the ACME Land Model (ALM) employs a non-unified treatment of hydrologic processes in the subsurface. Presently, ALM simulates transport of water in the subsurface via a theta-based Richards equation that is coupled to an unconfined aquifer model to account for groundwater-soil water interactions. This ad hoc treatment of vadose-phreatic zone interaction can result in unphysical model predictions of unsaturated soil layers below the predicted water table. A variably saturated flow model (VSFM) is required to overcome this shortcoming. We have developed a VSFM in ALM that uses the Portable, Extensible Toolkit for Scientific Computation (PETSc). Results of the VSFM model for various benchmark problems are presented.