Developing accurate, stable, and thermodynamically consistent numerical methods to simulate two-phase flows is critical for many applications. We develop numerical methods to solve thermodynamically consistent Cahn-Hilliard Navier-Stokes equations to simulate two-phase flows with deforming interfaces at various density contrasts. We develop three essentially unconditionally energy-stable time integration schemes. The first two time-integration schemes are fully implicit based on the pressure-stabilization technique. The third approach utilizes the projection method to decouple the pressure to improve the efficiency of the fully implicit scheme. We rigorously prove the energy stability of the time-discrete numerical schemes for the approaches with pressure-stabilization approaches. We also prove the existence of solutions of the advective-diffusive Cahn-Hilliard operator. We use a conforming continuous Galerkin (CG) finite element method in space equipped with a residual-based variational multiscale (RBVMS) procedure to stabilize the pressure in the first two approaches. In the third approach, we present a projection based framework extending the fully implicit method to a block iterative, hybrid semi-implicit-fully-implicit in time method. We use a semi-implicit time discretization for Navier-Stokes and a fully-implicit time discretization for Cahn-Hilliard equations. Pressure is decoupled using a projection step resulting in two linear positive semi-definite systems for velocity and pressure instead of the saddle point system of a pressure-stabilized method. All the resulting linear systems are solved using the efficient and scalable algebraic multigrid (AMG) method. We deploy all three approaches on a massively parallel numerical implementation using fast octree-based adaptive meshes in a computational framework called "Proteus". We perform a detailed scaling analysis of all three solvers. A comprehensive set of numerical experiments showing detailed comparisons with results from the literature for canonical cases are used to validate the methods for an extensive range of density ratios. This presentation is an overview of the main developments of Khanwale's PhD research. Khanwale will also discuss his plans to push the boundaries of phase-field methods and develop robust multiphysics solvers which couple scalar transport, electroconvection with multiphase flows.