The simulation of flow and transport in porous media such as aquifers often involve dealing with complex heterogeneities. They are characterized by varying hydrogeological properties which differ strongly from the adjacent medium and often lead to significant changes of the flow behavior. However detailed information about the location of such heterogeneities is not always known. The deterministic models thus need to be extended stochastically to quantify uncertainties. As mathematical model we use the capillarity‐free fractional flow formulation for two immiscible and incompressible fluid phases in a two‐dimensional and partitioned domain. To cope with the randomly located heterogeneity interfaces we employ a stochastic Galerkin (SG) method [4]. The physical space of this system then is modelled by a central upwind finite volume scheme [5] in combination with mixed finite elements [7]. (© 2016 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)