Abstract.  A novel numerical algorithm for computing incompressible, discontinuous, two-phase flows in two-dimensional, inhomogeneous, and isotropic porous media is presented. The algorithm uses Colella et al.’s hybrid sequential semi-implicit approach⁶ for both accuracy and efficiency of the calculations. The explicit part uses a high-order Godunov¹¹ scheme with a modified Van Leer geometrical slope limiter, similar to those used in shock dynamics. The implicit part is a two-step solver: the first step is a Crank-Nicolson saturation solver and the second one is a Poisson solver for the phase pressure. Both use fast multilevel multigrid solvers with the number of operations of the order of O[Nlog(N)], where N is the number of grid points. Two numerically stiff reservoir-engineering problems are presented to demonstrate the low numerical dispersion and second-order accuracy of our method.