We develop an accurate and robust numerical scheme for solving the incompressible hydrodynamically coupled Cahn–Hilliard system of the two-phase fluid flow system on complex surfaces. Our algorithm leverages a number of efficient techniques, including the subdivision-based isogeometric analysis (IGA) method for spatial discretization, the explicit Invariant Energy Quadratization (EIEQ) method for linearizing nonlinear potentials, the Zero-Energy-Contribution (ZEC) method for decoupling, and the projection method for the Navier–Stokes equation...