  Multiphase flow in porous media occurs in several disciplines including
petroleum reservoir engineering, petroleum systems' analysis, and CO$_2$
sequestration. While simulations often use a fully implicit discretization to
increase the time step size, restrictions on the time step often exist due to
non-convergence of the nonlinear solver (e.g. Newton's method). Here this
problem is addressed for the Buckley-Leverett equations, which model
incompressible, immiscible, two-phase flow with no capillary potential. The
equations are recast as a gradient flow using the phase-field method, and a
convex energy splitting scheme is applied to enable large timesteps, even for
high degrees of heterogeneity in permeability and viscosity. By using the
phase-field formulation as a homotopy map, the underlying hyperbolic flow
equations can be solved with large timesteps. For a heterogeneous test problem,
the new homotopy method allows the timestep to be increased by more than six
orders of magnitude relative to the unmodified equations while maintaining
convergence.
