Earth is characterized by episodes of glaciations and periods of minimal or no ice through geologic time. Using the linear energy balance model (EBM), nonlinear EBM with empirical ice sheet schemes, the general circulation model coupled with an ice sheet model, this study investigates the occurrence and stability of glaciations in geologic time. The simulations since the last glacial maximum (LGM) suggest that the summertime thawline of ice sheets conforms closely to the equatorward edge of the ice sheets and implies the relative stability toward deglaciation. CO2 levels are indispensable in controlling the initiation of ice sheet in the Cretaceous. At low CO2 levels, ice sheets exist in all periods no matter LGM or the last interglacial (LIG) orbital elements; however, at high CO2 levels ice sheets rarely exist. The simulations agree well with recent geological evidence of the hysteresis of glaciations in the Permo-Carboniferous. Gondwanaland reached its glacial maximum when CO2 level was roughly the same or slightly higher than the preindustrial value. With a further increase of CO2, deglaciation dominates and results in an ice free state. Again, if CO2 decreased to the present level, Gondwanaland would be glaciated once more and start a new cycle of glaciation and deglaciation. Simulations from five paleogeography maps in Gondwanaland with a suite of CO2 levels and different orbital elements reveal that paleogeography, CO2 levels and the Milankovitch cycles all contribute to the glaciations of Gondwanaland. This study shows that orbital elements alone are insufficient to account for the evolution of ice sheets. Net radiative forcing caused by greenhouse gases, such as CO2 and solar constant change are the primary drivers to glacial inception or demise. Continental geography, CO2 levels, solar constant change, and the Milankovitch cycles complicate the glacial history of Earth.