  We present a fully-coupled, implicit-in-time framework for solving a
thermodynamically-consistent Cahn-Hilliard Navier-Stokes system that models
two-phase flows. In this work, we extend the block iterative method presented
in Khanwale et al. [\textit{Simulating two-phase flows with thermodynamically
consistent energy stable Cahn-Hilliard Navier-Stokes equations on parallel
adaptive octree based meshes}, J. Comput. Phys. (2020)], to a fully-coupled,
provably second-order accurate scheme in time, while maintaining
energy-stability. The new method requires fewer matrix assemblies in each
Newton iteration resulting in faster solution time. The method is based on a
fully-implicit Crank-Nicolson scheme in time and a pressure stabilization for
an equal order Galerkin formulation. That is, 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. We deploy
this approach on a massively parallel numerical implementation using parallel
octree-based adaptive meshes. We present comprehensive numerical experiments
showing detailed comparisons with results from the literature for canonical
cases, including the single bubble rise, Rayleigh-Taylor instability, and
lid-driven cavity flow problems. We analyze in detail the scaling of our
numerical implementation.
