The speciation of calcium carbonate in water is important to the geochemistry of the world's oceans and has ignited significant debate regarding the mechanism by which nucleation occurs. Here, it is vital to be able to quantify the thermodynamics of ion pairing versus higher order association processes in order to distinguish between possible pathways. Given that it is experimentally challenging to quantify such species, here we determine the thermodynamics for ion pairing and multiple binding of calcium carbonate species using bias-enhanced molecular dynamics. In order to examine the uncertainties underlying these results, we derived a new polarizable force field for both calcium carbonate and bicarbonate in water based on the AMOEBA model to compare against our earlier rigid ion model, both of which are further benchmarked against ab initio molecular dynamics for the ion pair. Both force fields consistently indicate that the association of calcium carbonate ion pairs to form larger species is stable, though with an equilibrium constant that is lower than for ion pairing itself.