A two‐dimensional finite element model based on Galerkin's weighted residual approach and an upstream weighting technique was developed to predict simultaneous flow of water and oil in a three‐fluid phase system with gas assumed at constant pressure. Element matrices were computed using the influence coefficient method for both Picard and Newton‐Raphson nonlinear iteration schemes. A number of hypothetical simulations were performed in both one and two dimensions to evaluate the accuracy and efficiency of the various schemes with respect to handling of nonlinear soil properties, time marching, mass balance errors, soil nonhomogeneity, and effects of upstream weighting. Results indicate that the Picard scheme appears to be as effective as the Newton‐Raphson scheme while requiring substantially less computational effort if upstream weighting is employed. The Picard method without upstream weighting did not provide satisfactory convergence behavior. Results from problems involving extreme soil nonhomogeneity indicate that accurate solutions can be obtained with the Picard method with the proper use of upstream weighting. Saturation‐pressure derivative terms in the formulation can be evaluated using both analytical and chord‐slope methods. Results indicate that time averaging of these terms is critical for good mass balance results especially during redistribution. Overall mass balance is also sensitive to the fluid‐dependent scaling parameters of the constitutive