During recent years the authors and collaborators have been involved in an activity related to the constmction and analysis of large time step operator sphtting algorithms for the numerical simulation of multi-phase fiow in heterogeneous porous media. The purpose of these lecture notes is to review some of this activity. We illustrate the main ideas behind these novel operator sphtting algorithms for a basic two-phase flow model. Special focus is posed on the numerical solution algorithms for the saturation equation, which is a convection dominated, degenerate convection-diffusion equation. Both theory and applications are discussed. The gen eral background for the reservoir flow model is reviewed, and the main features of the numerical algorithms are presented. The basic mathematical results supporting the numerical algorithms are also given. In addition, we present some results from the BV solution theory for quasilinear degenerate parabolic equations, which provides the correct mathematical framework in which to analyse our nmnerical algorithms. Two-and three-dimensional munerical test cases are pre sented and discussed. The main conclusion drawn from the numerical experiments is that the operator sphtting algorithms indeed exhibit the property of resolving accurately internal layers with steep gradients, give very little numerical diffusion, and, at the same time, permit the use of large time steps. In addition, these algorithms seem to capture all potential combinations of convection and diffusion forces, ranging from convection dominated problems (including the pure hyperbolic case) to more diffusion dominated problems. CONTENTS 1.