Oxygen vacancies at ceria (CeO2) surfaces play an essential role in catalytic applications. However, during the past decade, the near-surface vacancy structures at CeO2(111) have been questioned due to the contradictory results from experiments and theoretical simulations. Whether surface vacancies agglomerate, and which is the most stable vacancy structure for varying vacancy concentration and temperature, are being heatedly debated. By combining density functional theory calculations and Monte Carlo simulations, we proposed a unified model to explain all conflicting experimental observations and theoretical results. We find a novel trimeric vacancy structure which is more stable than any other one previously reported, which perfectly reproduces the characteristics of the double linear surface oxygen vacancy clusters observed by STM. Monte Carlo simulations show that at low temperature and low vacancy concentrations, vacancies prefer subsurface sites with a local (2 × 2) ordering, whereas mostly linear surface vacancy clusters do form with increased temperature and degree of reduction. These results well explain the disputes about the stable vacancy structure and surface vacancy clustering at CeO2(111), and provide a foundation for the understanding of the redox and catalytic chemistry of metal oxides.