Numerical Solvers for the Richards Equation: A Comparative Review of Stability, Efficiency and Mass Conservation

Richards’ equation (RE) describes transient water flow in variably saturated porous media and plays a fundamental role in hydrological, agricultural, and environmental modelling. Its inherent nonlinearity, potential degeneracy, and sensitivity to hydraulic parameterization make accurate numerical solutions a persistent challenge. This review offers a systematic comparative analysis of numerical methods for the RE across three key dimensions: spatial discretization, primary variable formulation, and iterative linearization strategy. Finite Difference, Finite Element, and Finite Volume methods are evaluated for their theoretical properties and performance under demanding conditions such as sharp infiltration fronts and heterogeneous parameter fields. Meshless methods are also assessed for their geometric flexibility and accuracy. Iterative linearization schemes are compared in terms of convergence, robustness, and computational cost. The analysis is structured around three performance metrics: numerical stability, computational efficiency, and mass conservation. Mass conservation is found to depend critically on both the primary variable formulation and the discretization framework, with mixed and finite volume approaches offering the strongest guarantees. The review concludes that no universally optimal solver exists; method selection must be tailored to specific soil conditions, boundary conditions, and computational constraints. Key open challenges include convergence under extremely dry conditions, multiscale heterogeneity, and coupled flow-transport processes.