To study the dependence of the relative stability of surface (V A ) 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 the one in the surface stability of the facets, i.e. (111) < (110) < (100). The results also reveal 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 V A is more stable than the V B at the (110) facet whereas at the (111), it is the other way around, and at the (100) facet, both the V A and the VB have similar stability. The stability of the V A vacancies, compared to that of the V B , 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, as well as the repulsive Ce3+−Ce3+ interactions.