Flow in the subsurface, which is relevant in such fields as geohydrology and petroleum engineering, deals with complex physical processes, e.g., advection, capillarity and gravity. The proper description of these processes and their relative importance for the output of interest are challenging tasks which need to be addressed if one is to achieve accurate and fast simulations of flow in such porous media. In addition, geological heterogeneities further increase the complexity of the problems and unstructured grids are often needed to characterize complex natural formations. So one of the main objectives of current research is to develop robust numerical algorithms which reduce the computational burden of realistic, large-scale problems and at the same time ensure reliable solutions. This is also the goal of this dissertation. Here, the main focus is on the fractional flow formulation of two-phase flow in porous media, in which robust discretization methods for the pressure and the saturation equation are studied separately.For the pressure equation, the multi-point flux approximation (MPFA) L-method is investigated including the Dirichlet boundary influence, numerical studies, geometrical interpretations of the L-method for homogeneous porous media and optimal H 1 and L 2 error estimates for the cases with homogeneous media on uniform grids. The second main part of the work is the investigation and development of time and space discretizations for the saturation equation based on the operator splitting (OS) concept in combination with the streamline approach. Special focus is put on the corrected operator splitting (GOS) method which allows for large time steps. The novel aspects of this part are the developments of a two-scale OS approach which solves different physical processes on different grids and a new GOS method which provides a promising tool for efficiently simulating complex porous media flow with gravity.