To study the dependence of the relative stability of surface (VA) and subsurface (VB) oxygen vacancies with the crystal facet of CeO2, the reduced (100), (110) and (111) surfaces, with two different concentrations of vacancies, were investigated by means of density functional theory (DFT+U) calculations. The results show that the trend in the near-surface vacancy formation energies for comparable vacancy spacings, i.e. (110) < (100) < (111), does not follow that in the surface stability of the facets, i.e. (111) < (110) < (100). The results also revealed that the preference of vacancies for surface or subsurface sites, as well as the preferred location of the associated Ce3+ polarons, are facet and concentration dependent. At the higher vacancy concentration, the VA is more stable than the VB at the (110) facet whereas at the (111), it is the other way around, and at the (100) facet, both the VA and the VB have similar stability. The stability of the VA vacancies, compared to that of the VB, is accentuated as the concentration decreases. Nearest neighbor polarons to the vacant sites are only observed for the less densely packed (110) and (100) facets. These findings are rationalized in terms of the packing density of the facets, the lattice relaxation effects induced by vacancy formation and the localization of the excess charge, and the repulsive Ce3+ - Ce3+ interactions.