Density of seawater is a critical property that controls ocean dynamics. Previous works suggest the use of the δ18O calcite of foraminifera as a potential proxy for paleodensity. However, potential quantitative reconstructions were limited to the tropical and subtropical surface ocean and without an explicit estimate of the uncertainty in calibration model parameters. We developed the use of the δ18Oc of planktonic foraminifera as a surface paleodensity proxy using Bayesian regression models calibrated to annual surface density. Predictive performance of the models improves when we account for inter-species specific differences. We investigate the additional uncertainties that could be introduced by potential evolution of the δ18Oc-density relationship with time – from the last glacial maximum (LGM) to the preindustrial (PI) – through the combination of past isotope enabled climate model simulations and a foraminiferal growth module. We demonstrate that additional uncertainties are weak globally, except for the Nordic Seas region. We applied our Bayesian regression model to LGM and Late Holocene (LH) δ18Oc foraminifera databases to reconstruct annual surface density during these periods. We observe stronger LGM density value changes at low latitudes compared to mid latitudes. These results will be used to evaluate numerical climate models in their ability to simulate ocean surface density during the extreme climatic period of the LGM. The new calibration has great potential to reconstruct the past temporal evolution of ocean surface density over the Quaternary. Under climates outside the Quaternary period and in ocean basins characterized by anti-estuary circulation, like the current Mediterranean Sea and Red Sea, our calibration could provide density estimates with larger uncertainty, a point that requires further investigations.