Using a classical Newton-Krylov method to solve the resulting nonlinear system of two-phase flows in porous media often suffers from slow convergence or failure in line search. We propose two nonlinear elimination preconditioning strategies to handle this issue by performing subspace correction to remove the local strong nonlinearities. Numerical experiments show that the proposed methods are more robust and faster than the existing method with respect to some physical and numerical parameters, and scalable to thousands of processes.