Two-phase flow in porous media is modeled at large spatial scales by two coupled degenerate, nonlinear parabolic equations, or sometimes a single degenerate parabolic equation. These models arise from averaging microscale free-boundary problems, with both modeling approaches relevant for many natural and industrial processes. We consider specifically water flow in the shallow and variably saturated subsurface environments, such as levees, storm water infiltration systems, and natural coastal environments. Due to the degeneracy of the equations and the density and viscosity ratios for air/water flows, robust, accurate, and efficient discretizations on unstructured meshes have been difficult to extend to practical higher-order discretizations. Recent advances in limiting of higher order continuous finite element methods, alongside advances in analysis of degenerate equations, has opened up routes to genuine advances in unsaturated flow modeling. This work presents work toward higher-order local polynomial approximations for air/water flow, including approaches based on the Richards' equation extending bound-preserving low-order schemes to higher-order methods.