Summary: In this paper we present PDE and finite element analyses for a system of PDEs consisting of the Darcy equation and the Cahn–Hilliard equation, which arises as a diffuse interface model for the two-phase Hele-Shaw flow. In the model the two sets of equations are coupled through an extra phase induced force term in the Darcy equations and a fluid induced transport term in the Cahn–Hilliard equation. We propose a fully discrete implicit finite element method for approximating the PDE system, which consists of the implicit Euler method combined with a convex splitting energy strategy for the temporal discretization, the standard finite element discretization for the pressure, and a split (or mixed) finite element discretization for the fourth-order Cahn–Hilliard equation. It is shown that the proposed numerical method satisfies a mass conservation law in addition to a discrete energy law that mimics the basic energy law for the Darcy–Cahn–Hilliard phase field model and holds uniformly in the phase field parameter \(\varepsilon\). With the help of the discrete energy law, we first prove that the fully discrete finite method is unconditionally energy stable and uniquely solvable at each time step. We then show that, using the compactness method, the finite element solution has an accumulation point that is a weak solution of the PDE system. As a result, the convergence result also provides a constructive proof of the existence of global-in-time weak solutions to the Darcy–Cahn–Hilliard phase field model in both two and three dimensions. Numerical experiments based on the overall solution method of combining the proposed finite element discretization and a nonlinear multigrid solver are presented to validate the theoretical results and to show the effectiveness of the proposed fully discrete finite element method for approximating the Darcy–Cahn–Hilliard phase field model.