The dynamical ellipticity of a planet expresses the departure of its mass distribution from spherical symmetry. It enters as a parameter in the description of a planet's precession and nutation, as well as other rotational normal modes. In the case of the Earth, uncertainties in this quantity's history produce an uncertainty in the solutions for the past evolution of the Earth-Moon system. Constraining this history has been a target of interdisciplinary efforts as it represents an astro-geodetic parameter whose variation is shaped by geophysical processes, and whose imprints can be found in the geological signal. We revisit the classical problem of its variation during ice ages, where glacial cycles exerted a varying surface loading that had altered the shape of the geoid. In the framework of glacial isostatic adjustment, and with the help of a recent paleoclimatic proxy of ice volume, we present the evolution of the dynamical ellipticity over the Cenozoic ice ages. We map out the problem in full generality identifying major sensitivities to surface loading and internal variations in parameter space. This constrained evolution is aimed to be used in future astronomical computations of the orbital and insolation quantities of the Earth.