The spatial extent of the Sahara (the largest nonpolar desert) has significant impacts on the livelihood of people residing in its surrounding areas. Despite the fact that climate change would foreseeably impact the location and size of the desert, its future responses (i.e., advance or retreat) are rarely explored in previous studies. Here, through the development of an ensemble Bayesian discriminant analysis approach, we use 10 of the latest high-resolution GCM (global climate model) simulations to document robust annual and seasonal responses of the Sahara Desert to twenty-first-century climate change, with the consideration of modeling uncertainties. We find northward shifts of the Sahara/Sahel and eastern expansion of the nondesert zone under both SSP2–4.5 and SSP5–8.5 scenarios, the former more pronounced in the wet season and the latter in the dry season. Countries located near the Mediterranean may thus experience higher risks of drought, while the projected retreat of the Saharan southern boundary will be beneficial to the local water availability of proximal countries. Significance Statement Given that sub-Saharan Africa is one of the most vulnerable regions to climate change, the Sahara’s expansion would bring unexpected health risks to billions of people. It is thus vital to understand its robust response to global warming. However, previous studies are merely focused on using a simple precipitation threshold as the definition criterion to estimate the varying size of the Sahara Desert. In addition, significant uncertainty in precipitation projections also limits relevant investigations of the Sahara’s future responses. Here, by developing an ensemble Bayesian discriminant analysis approach, we could provide an objective basis for desert identification under large intermodel uncertainty. Further, we find significant northward shifts of both the Sahara and the Sahel, which may induce higher risks of drought over the northwest of North Africa.