We develop a novel entropy–stable discontinuous Galerkin approximation of the incompressible Navier–Stokes/Cahn–Hilliard system for p–non–conforming elements. This work constitutes an evolution of the work presented by Manzanero et al. ((2020) [10]), as it extends the discrete analysis into supporting p–adaptation (p–refinement/coarsening). The scheme is based on the summation–by–parts simultaneous–approximation term property along with Gauss–Lobatto points and suitable numerical fluxes. The p–non–conforming elements are connected through the classic mortar method, the use of central fluxes for the inviscid terms, and the BR1 scheme with additional dissipation for the viscous fluxes. The scheme is proven to retain its properties of the original conforming scheme when transitioning to p–non–conforming elements and to mimic the continuous entropy analysis of the model. We focus on dynamic polynomial adaptation as the applications of interest are unsteady multiphase flows. In this work, we introduce a heuristic adaptation criterion that depends on the location of the interface between the different phases and utilises the convection velocity to predict the movement of the interface. The scheme is verified to be total phase conserving, entropy–stable and freestream preserving for curvilinear p–non–conforming meshes. We also present the results for a rising bubble simulation and we show that for the same accuracy we get a ×2 to ×6 reduction in the degrees of freedom and a 41% reduction in the computational time. We compare our results for the three–dimensional dam break test case against experimental and numerical data and we show that a ×4.3 to ×9.5 reduction of the degrees of freedom and a 51% reduction in the computational time can be achieved compared to the p–uniform solution.