  In this work we present a mass conservative numerical scheme for two-phase
flow in porous media. The model for flow consists on two fully coupled,
non-linear equations: a degenerate parabolic equation and an elliptic equation.
The proposed numerical scheme is based on backward Euler for the temporal
discretization and mixed finite element method (MFEM) for the discretization in
space. Continuous, semi-discrete (continuous in space) and fully discrete
variational formulations are set up and the existence and uniqueness of
solutions is discussed. Error estimates are presented to prove the convergence
of the scheme. The non-linear systems within each time step are solved by a
robust linearization method. This iterative method does not involve any
regularization step. The convergence of the linearization scheme is rigorously
proved under the assumption of a Lipschitz continuous saturation. The case of a
H\"older continuous saturation is also discussed, a rigorous convergence proof
being given for Richards' equation. Numerical results are presented to sustain
the theoretical findings.
