  The paper considers a thermodynamically consistent phase-field model of a
two-phase flow of incompressible viscous fluids. The model allows for a
non-linear dependence of fluid density on the phase-field order parameter.
Driven by applications in biomembrane studies, the model is written for
tangential flows of fluids constrained to a surface and consists of (surface)
Navier-Stokes-Cahn-Hilliard type equations. We apply an unfitted finite element
method to discretize the system and introduce a fully discrete time-stepping
scheme with the following properties: (i) the scheme decouples the fluid and
phase-field equation solvers at each time step, (ii) the resulting two
algebraic systems are linear, and (iii) the numerical solution satisfies the
same stability bound as the solution of the original system under some
restrictions on the discretization parameters. Numerical examples are provided
to demonstrate the stability, accuracy, and overall efficiency of the approach.
Our computational study of several two-phase surface flows reveals some
interesting dependencies of flow statistics on the geometry.
