The scale of heterogeneities in reservoirs is often smaller than the grid size used in large-scale reservoir simulations. Relative permeabilities have the foremost effect on fluid flow and small-scale fractional flow must be upscaled to the grid size in flow simulations. Thus, systems with fractional flow heterogeneities, i.e. relative permeability variations, have to be solved and effective relative permeability functions determined. In this paper, benchmark analytical solutions are developed for a system consisting of two media in series where each medium is characterized with a different set of relative permeabilities, residual saturations and porosities. The analytical solutions show a significant discontinuity in the saturation and concentration profiles at the interface of the two media. Numerical results using several weighting schemes are compared against the analytical solutions for two-phase immiscible and partially miscible systems. It is shown that single-point upstream weighting requires about 100 grid blocks to capture the discontinuity at the interface, whereas third-order TVD weighting requires much fewer. Lastly, the validity of the JBN method and harmonic averaging for determination of effective relative permeabilities and overall pressure drop is tested. Both the JBN method and harmonic averaging cannot reproduce the pressure drop across the composite media prior to water breakthrough.