The paper investigates the numerical behavior of standard approximation methods for the unsaturated flow equation. Solution using the h-based formulation and a backward Euler time discretization is shown to produce unacceptably large mass balance errors for many example calculations. This is true for any iteration method (Picard, Newton-Raphson, etc.). It is also true for both finite difference and finite element approximations in space, although finite elements are generally inferior to finite differences. A modified numerical approach is proposed that alleviates the mass balance problems discussed above. This approach is based on a fully implicit (backward Euler) time approximation applied to the mixed form of the unsaturated flow equation.