The anomalous volume isotope effect (VIE) of ice Ih is calculated and analyzed based on the quasi-harmonic approximation to account for nuclear quantum effects in the Helmholtz free energy. While a lot of recently developed polarizable many-body potential functions give a normal VIE contrary to experimental results, we find that one of them, MB-pol, yields the anomalous VIE in good agreement with the most recent high-resolution neutron diffraction measurements─better than DFT calculations. The short-range three-body terms in the MB-pol function, which are fitted to CCSD(T) calculations, are found to have a surprisingly large influence. A vibrational mode group decomposition of the zero-point pressure together with a hitherto unconsidered benchmark value for the intramolecular stretching modes of H2O ice Ih obtained from Raman spectroscopy data unveils the reason for the VIE: a delicate competition between the latter and the librations.