Abstract The Fully Implicit Method (FIM) is widely used to discretize the nonlinear conservation equations that govern multiphase flow in porous media. Although FIM is unconditionally stable, the nonlinear solver, which is typically based on Newton’s method, often fails to converge for large time steps. As the time step size increases, the complexity of the nonlinear coupling increases. In the presence of strong buoyancy in a heterogeneous domain, for example, counter-current flow may be present, and phase-flow reversal - as a function of time or Newton iteration - becomes an issue. It does not take too many interfaces with flow reversal across them to make the nonlinear solver stall. Even when the Newton-based iterative procedure converges, the cost associated with updating all the unknowns simultaneously - everywhere and all the time - can be quite expensive. Conventional sequential-implicit strategies can be used to reduce the cost, but they suffer from severe restrictions on the allowable timestep size. We propose a new nonlinear solution strategy to improve the robustness and efficiency of FIM. Specifically, for each FIM Newton iteration, we update pressure by solving a linear system obtained from the full Jacobian. Then, we construct a directed graph with grid-cells as nodes and phase fluxes at grid-cell interfaces as edges, to identify the counter-current flow regions. After that, we employ a topological ordering that makes it possible to visit the components once, from upstream to downstream. The graph components vary in size from a single cell to multiple ‘strongly connected’ cells. For each component made up of a single cell, we update the saturation nonlinearly. For each component made up of multiple cells (strongly connected region), we perform a simultaneous linear update of both saturation and pressure. The saturation iterates are safeguarded based on ‘trust regions’ of the flux (fractional flow) functions. We compare the new solver with the current state-of-the-art for a wide range of heterogeneous problems with special focus on counter-current flow due to buoyancy. Our results show that the proposed FIM nonlinear solver converges for large time steps, with excellent convergence rates and lower computational cost per iteration compared with standard or safeguarded Newton.