In this article, we study a phase field model for a two-layer fluid where the temperature dependence of both the density (buoyancy forces) and the surface tension (Marangoni effects) is considered. The phase field model consisting of a modified Navier–Stokes equation, a Cahn–Hilliard phase field equation and an energy transport equation is derived through an energetic variational procedure. An appropriate variational form and a continuous finite element method are adopted to maintain the underlying energy law to its greatest extent. A few examples for Bénard–Marangoni convection in an Acetonitrile and n-Hexane two-layer fluid system heated from above will be computed to justify our phase field model and further show the good performance of our methods. In addition, an interesting experiment will be performed to show the competition between the Marangoni effects and the buoyancy forces.