The extension of Darcy's law to two-phase flows is central to modeling several critical processes in civil, environmental, and coastal engineering applications, including infiltration and seepage. These processes are inherently multiscale, resulting in highly nonlinear closure relations and material heterogeneity at engineering scales. Early work extending higher-order methods to realistic model equations encountered many challenges, including avoiding mass conservation errors in time and recovering suitable local mass conservation properties for numerical velocity fields in space. Through several decades, numerical models for these processes have remained notorious for their computational expense and lack of robust solvers. Recent work reviving Flux-Corrected Transport (FCT) and related methods to build numerical schemes that preserve fundamental physical properties has interesting extensions to multiphase flow models at engineering scales. In this work we consider the extension of FCT-type schemes to VMS schemes for unstructured finite element methods to preserve local saturation bounds and maintain mass conservation in time and space, while still achieving optimal order. By avoiding low-order regularizations, such as nonlinear shock-capturing, the schemes maintain higher-order convergence and reduce the number of tunable parameters facing end users. The algebraic stabilization approaches inherited from FCT enable these properties and provide other useful properties as well, such as robustness for highly anisotropic meshes and materials. We consider convergence results on a range of benchmark problems as well as application to green infrastructure.