For many reservoir simulation problems, critical flow gradients occupy only a small percentage of the flow domain, and the ability to perform local mesh refinement is highly desirable. Moving away from conventional grid systems raises many issues in terms of flow equation discretization. Hyperbolic and elliptic adaptive-grid-interface flux approximations are presented and discretization issues are discussed. Aquifer remediation simulations involving transport of steep saturation gradients demonstrate that benefits of high resolution and reduced computational cost are achievable when coupling adaptive local grid refinement with higher-order schemes in three dimensions.