Oxygen vacancies on ceria (CeO₂) surfaces play a crucial role in catalytic applications, yet whether vacancies are at surface or subsurface sites on reduced CeO₂(111), and whether vacancies agglomerate or repel each other, is still under discussion, with few and inconsistent experimental results. By combining density-functional theory (DFT) in the DFT + U (U is an effective onsite Coulomb interaction parameter) approach and statistical thermodynamics, we show that the energetically most stable near-surface oxygen vacancy structures for a broad range of vacancy concentrations, θ (1/16 ≥ θ ≥ 1 monolayer) have all vacancies at subsurface oxygen sites and predict that the thermodynamically stable phase for a wide range of reducing conditions is a (2 × 2) ordered subsurface vacancy structure (θ = ¼). Vacancy-induced lattice relaxations effects are crucial for the interpretation of the repulsive interactions, which are at the basis of the vacancy spacing in the (2 × 2) structure. The findings provide theoretical data to support the interpretation of the most recent experiments, bringing us closer to solving the debate.