Implementation and Evaluation of Fully 3D Richards' Equation in MODFLOW 6

This paper presents the development and benchmarking of a fully three-dimensional (3D) Richards' equation solver integrated into MODFLOW 6. The popular unsaturated zone package (the UZF package) currently in MODFLOW relies on simplified approaches that neglect capillary effects and assumes one-dimensional (1D) vertical flow. In contrast, the Richards' equation implementation presented here is based on a 3D head-based mixed formulation using the van Genuchten moisture retention and relative permeability functions to capture the more detailed physics of variably saturated flow. The numerical implementation within the MODFLOW 6 framework leverages the advanced, robust numerical techniques already built into the software including upstream weighting, Newton-Raphson and Picard linearization approaches, advanced linear solvers and a flexible finite-volume framework. Three benchmark tests were used to evaluate the performance of this Richards' equation implementation in MODFLOW 6 against established codes such as HYDRUS, MODFLOW VSF, SUTRA, and USG-T. Results indicate that this package accurately reproduces moisture dynamics and capillary pressure distributions with computational efficiency and robustness. This advancement not only improves the physical realism of saturated/unsaturated subsurface flow simulations but also provides a practical tool for further research in variably saturated flow modeling within the MODFLOW framework.