A variational inequality formulation of two-phase flow problem is proposed. The variational inequality ensures the physical feasibility saturation. A fully implicit method with adaptive dt is applied to achieve high stability. A nonlinear elimination preconditioner is used to accelerate nonlinear solve. Numerical results show the proposed approach is more robust than classical IMPES.