The U.S. Geological Survey ground-water flow model MODFLOW often has difficulties solving freesurface conditions due to the highly non-linear governing equations. MODFLOW currently handles freesurface conditions by “wetting” and “drying” cells. A cell becomes “dry” or inactive when the simulated water level falls below the base elevation of the cell. Likewise, a cell is allowed to “wet” or become activewhen water levels in adjacent cells are above a user-defined threshold. This approach can be numerically unstable, especially in the presence of heterogeneity. A new algorithm is developed which allows flow along diagonal grid points near the free-surface boundary. Geometric relations are used to estimate the average hydraulic conductivity for diagonal cells. The results indicate that the discontinuity that typically develops when the water table crosses a cell interface is effectively removed when flow is allowed along diagonal water table cells. The diagonal-cell approach is also stable under heterogeneous conditions. The drawback to the approach is that it imposes additional bands in the coefficient matrix that is used to solve the linear system. Because MODFLOW solvers are designed to handle heptadiagonal (seven bands) coefficient matrices, solver modifications would be required for this method to be implemented into MODFLOW.