So-called corrected operator splitting methods are applied to a 1-D scalar advection diffusion equation of Buckley-Leverett type with general initial data. Front tracking and a 2nd order Godunov method are used to advance the solution in time. Diffusion is modelled by piecewise linear finite elements at each new time level. To obtain correct structure of shock fronts independently of the size of the time step, a dynamically defined residual flux term is grouped with diffusion. Different test problems are considered, and the methods are compared with respect to accuracy and runtime. Finally, wc extend the corrected operator splitting to 2-D equations by means of dimensional splitting, and wc apply it to a Buckley-Leverett problem including gravitational effects. KARLSEN, BRUSDAL, DAHLE, EVJE, AND LIE simple analytic expressionThe capillary diffuse coefficient u(u) is generally a nonlinear (bell-shaped) funct lOn of u, which we recreate using the expression We note that the diffusion coefficient becomes zero at u = 0, 1, so that (1) is an example of a degenerate parabolic equation; see Wpert and Hudjaev [37] for global existence, uniqueness, and stability results in BV space for the initial value problem with general "(") > 0. and Zhuo-qun and Jun-yu [42] for similar results concerning the initial-boundary value problem. Properties of BV solutions, such as regularity, have been studied by Jun-Ning [25]. For further results on degenerate equations we refer to a recent survey paper by Chen [6] and the references therein.The scaling parameter e in front of the capillary diffusion term is usually small for reasonable flow rates. Consequently, this term can be neglected if the main problem is to trace fluid interfaces in which case the problem is reduced to solving a hyperbolk conservation law. However, in many applications some detailed information on the structure of fronts is important, and the diffusion term cannot be neglected. Since we believe that it is important to obtain correct placement and