Groundwater contamination previously occurred at a broad range of locations in present-day China. There are thousands of kinds of contaminants which can be divided into soluble and insoluble categories in groundwater. In recent years, the non-aqueous phase liquid (NAPL) pollution that belongs to the multi-phase seepage flow phenomenon has become an increasingly prominent topic due to the challenge brought by groundwater purification and its treatment. Migrating with seepage flow and moving into the potable water sources, these contaminants directly endanger people’s health. Therefore, it is necessary to research how these contaminants not only migrate, but also are then accordingly remedied. First, as an analysis means, an effective numerical method is necessary to be built. A three-dimensional finite element method program for analyzing two-phase flow in porous media, which can be applied to the immiscible contaminant transport problem in subsurface flow has been developed in this paper. The fundamental theory and numerical discretization formulations are elaborated. The numerical difficulty brought about by the distinct non-linearity of the temporal evolution of saturation-dependent variables is overcome by the mixed-form formulation. The effectiveness of simultaneous solution (SS) method and its improvement in efficiency are explained. Finally, two computational examples are given for verifying the correctness and demonstrating the preliminary applicability. In addition, the function of two-phase immiscible flow, especially in Fast Lagrangian Analysis of Continua (FLAC) is used to simulate the same examples and the results are compared to further verify the correctness of the numerical development.