A review of numerical methods for flow and transport modeling in the vadose zone

Accurate modeling of subsurface flow and transport processes in the vadose zone, governed by the Richards equation and advection‐dispersion equation, respectively, poses major numerical challenges due to the strong nonlinearity of the governing partial differential equations, spatial heterogeneity and anisotropy of the media parameters, and complex coupling between the physical processes. This paper presents an extensive review of the state‐of‐the‐art numerical modeling techniques, evaluating both traditional Eulerian frameworks and Lagrangian approaches. The paper quantitatively examines the global and local mass conservation characteristics of traditional numerical methods—finite difference, finite element, and finite volume (FV), through a comparative mass balance error analysis. The results demonstrate the high mass balance efficiency of the FV method, highlighting its suitability and robustness for highly nonlinear vadose zone modeling. Furthermore, this review work explores advanced modeling paradigms for capturing non‐Fickian transport phenomena resulting from the strong heterogeneity of natural porous media. The paper also critically assesses the emergence of artificial intelligence and machine learning (AI/ML) surrogates, specifically physics‐informed neural networks and operator learning, in hydrogeological modeling. Overall, this paper identifies critical future research avenues, emphasizing the necessity for fully coupled multi‐physics models (flow, mass, and heat transport) and the establishment of novel benchmark problems to simulate emerging contaminants, reactive transport, and complex hydrodynamics in data‐scarce subsurface environments.