Bound-Preserving Flux-Corrected Transport Methods for Solving Richards' Equation

Simulating infiltration in porous media using Richards' equation remains computationally challenging due to its parabolic structure and nonlinear coefficients. While a wide range of numerical methods for differential equations have been applied over the past several decades, basic higher-order numerical methods often fail to preserve physical bounds on water pressure and saturation, leading to spurious oscillations and poor iterative solver convergence. Instead, low-order, bound-preserving methods have been preferred. The combination of mass lumping and relative permeability upwinding preserves bounds but degrades accuracy to first order in space. Flux-corrected transport is a high-resolution numerical technique designed for combining the bound-preserving property of low-order schemes with the accuracy of high-order methods, by blending the two methods through limited anti-diffusive fluxes. In this work, we extend flux-corrected transport schemes to the nonlinear, degenerate parabolic structure of Richards' equation, verify attainment of second-order convergence on unstructured meshes, and demonstrate applications to stormwater management infrastructure.