  In this article we consider the numerical modeling and simulation via the
phase field approach of two-phase flows of different densities and viscosities
in superposed fluid and porous layers. The model consists of the
Cahn-Hilliard-Navier-Stokes equations in the free flow region and the
Cahn-Hilliard-Darcy equations in porous media that are coupled by seven domain
interface boundary conditions. We show that the coupled model satisfies an
energy law. Based on the ideas of pressure stabilization and artificial
compressibility, we propose an unconditionally stable time stepping method that
decouples the computation of the phase field variable, the velocity and
pressure of free flow, the velocity and pressure of porous media, hence
significantly reduces the computational cost. The energy stability of the
scheme effected with the finite element spatial discretization is rigorously
established. We verify numerically that our schemes are convergent and
energy-law preserving. Ample numerical experiments are performed to illustrate
the features of two-phase flows in the coupled free flow and porous media
setting.
