We present an anisotropic mesh adaptation procedure based on Riemannian metrics for the simulation of two-phase incompressible flows with non-matching densities. The system dynamics are governed by the Cahn-Hilliard Navier-Stokes (CHNS) equations, discretized with mixed finite elements and implicit time-stepping. The volume-averaged formulation of Abels et al. is considered [1]. Spatial accuracy is controlled throughout the simulation by the global transient fixed-point method from Alauzet et al. [2], in which the simulation time is divided into sub-intervals, each associated with an adapted anisotropic mesh. The simulation is run in a fixed-point loop until convergence of each mesh–solution pair. Each iteration takes advantage of the previously computed solution and accurately predicts the flow variations, so that the adapted mesh never lags behind the interface. This ensures that the mesh always captures the fluid-fluid interface, and allows for a dynamic control of the interface thickness at a fraction of the computational cost compared to uniform or isotropic grids. Moreover, using a reduced number of time sub-intervals reduces the transfer error from one mesh to another, which would otherwise eventually spoil the numerical solution. The overall adaptive procedure is verified with manufactured solutions and on well-known academic benchmarks, namely the rising bubbles and Rayleigh-Taylor instability.