Dewatered or "dry" grid cells in the USGS ground water modeling software MODFLOW may cause nonphysical artifacts, trigger convergence failures, or interfere with parameter estimation. These difficulties can be avoided in two dimensions by modifying the spatial differencing scheme and the iterative procedure used to resolve nonlinearities. Specifically, the spatial differencing scheme is modified to use the water level on the upstream side of a pair of adjacent cells to calculate the saturated thickness and hence intercell conductance for the pair. This makes it possible to explicitly constrain the water level in a cell to be at or above the cell bottom elevation without introducing nonphysical artifacts. Thus constrained, all initially active cells will remain active throughout the simulation. It was necessary to replace MODFLOW's Picard iteration method with the Newton-Raphson method to achieve convergence in demanding applications involving many dry cells. Tests using a MODFLOW variant based on the new method produced results nearly identical to conventional MODFLOW in situations where conventional MODFLOW converges. The new method is extremely robust and converged in scenarios where conventional MODFLOW failed to converge, such as when almost all cells dewatered. An example application to the Edwards Aquifer in south-central Texas further demonstrates the utility of the new method.