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. Recent studies have identified regions within a watershed that have a water table between 1 and 5 m below the ground surface to have the largest influence of groundwater dynamics on surface energy budgets. Numerous numerical studies have shown impacts of near-surface soil moisture dynamics on several key Earth system processes, including runoff, surface energy partitioning, vegetation dynamics, soil biogeochemistry, and greenhouse gas emissions. Despite the obvious need to accurately represent soil moisture dynamics, the current version of the Community Land Model (CLM) employs a non-unified treatment of hydrologic processes in the subsurface. To overcome this shortcoming, we implemented a variably saturated Richards equation (RE) that uses a water equation of state. The variably saturated RE differential equation is first rewritten as a system of differential algebraic equations (DAEs) and the Method of Lines (MOL) is then used to spatially discretize the system of DAEs. The resulting discretized DAEs system is evolved in time using a variable-order, variable-step-size backward difference formula. The Portable, Extensible Toolkit for Scientific Computation (PETSc) is used for numerically integrating the DAE system. Numerical results and comparisons with analytical solutions are presented for several test problems.