A two-phase Newtonian surface fluid is modelled as a surface Cahn–Hilliard–Navier–Stokes equation using a stream function formulation. This allows one to circumvent the subtleties in describing vectorial second-order partial differential equations on curved surfaces and allows for an efficient numerical treatment using parametric finite elements. The approach is validated for various test cases, including a vortex-trapping surface demonstrating the strong interplay of the surface morphology and the flow. Finally the approach is applied to a Rayleigh–Taylor instability and coarsening scenarios on various surfaces.