In reservoir simulation, solution of the coupled systems of nonlinear algebraic equations that are associated with fully-implicit (backward Euler) discretization is challenging. Having a robust and efficient nonlinear solver is necessary in order for reservoir simulation to serve as the primary tool for managing the recovery processes of large-scale reservoirs. However, there are several outstanding challenges that are intimately connected to the nature of the problem. Given a set of sources and sinks, the variation in the total velocity can span many orders of magnitude due to the extreme levels of permeability heterogeneity in large-scale subsurface porous formations. Moreover, multiple and complex saturation fronts need to be properly resolved throughout the three-dimensional reservoir model of interest. Add to that the fact that the numerical simulations need to predict the field-scale recovery performance over many decades, and the challenge of developing robust and efficient nonlinear solvers across a very wide parameter space becomes quite clear. Here, we develop a continuation method based on the use of a dissipation operator to simulate nonlinear two-phase flow and transport in heterogeneous porous media in the presence of viscous, gravitational, and capillary forces. The homotopy is constructed by adding numerical dissipation to the discrete conservation equations. A continuation parameter is introduced to control the amount of dissipation. Numerical evidence of multi-dimensional models and detailed analysis of single-cell problems are used to explain how the dissipation operator improves the nonlinear convergence of the coupled system of equations. An adaptive strategy to determine the optimum dissipation coefficient is proposed. The adaptive dissipation is computed locally for each cell interface, and is applicable to coupled flow and transport. We demonstrate the efficacy of the new nonlinear solver using several examples, including 1D scalar transport and 2D heterogeneous problems with fully-coupled flow and transport. The new solver exhibits superior convergence properties compared with the standard Newton solver used in reservoir simulation. We show that the new solver works robustly for a wide range of flow conditions without parameter tuning of the dissipation operator.