We present a phase-field model for multiphase flow for an arbitrary number of immiscible incompressible fluids with variable densities and viscosities. The model consists of a system of the Navier−Stokes equations coupled to multicomponent Cahn−Hilliard variational inequalities. The proposed formulation admits a natural energy law, preserves physically meaningful constraints and allows for a straightforward modelling of surface tension effects. We propose a practical fully discrete finite element approximation of the model which preserves the energy law and the associated physical constraints. In the case of matched densities we prove convergence of the numerical scheme towards a weak solution of the continuous model. The convergence of the numerical approximations also implies the existence of weak solutions. Furthermore, we propose a convergent iterative fixed-point algorithm for the solution of the discrete nonlinear system of equations and present several computational studies of the proposed model.