  We report on simulations of two-phase flows with deforming interfaces at
various density contrasts by solving thermodynamically consistent Cahn-Hilliard
Navier-Stokes equations. An (essentially) unconditionally energy-stable
Crank-Nicolson-type time integration scheme is used. Detailed proofs of energy
stability of the semi-discrete scheme and for the existence of solutions of the
advective-diffusive Cahn-Hilliard operator are provided. In space we discretize
with a conforming continuous Galerkin finite element method in conjunction with
a residual-based variational multi-scale (VMS) approach in order to provide
pressure stabilization. We deploy this approach on a massively parallel
numerical implementation using fast octree-based adaptive meshes. A detailed
scaling analysis of the solver is presented. Numerical experiments showing
convergence and validation with experimental results from the literature are
presented for a large range of density ratios.
