  We present a comparison of galaxy atomic and molecular gas properties in
three recent cosmological hydrodynamic simulations, Simba, EAGLE, and
Illustris-TNG, versus observations from $z\sim 0-2$. These simulations all rely
on similar sub-resolution prescriptions to model cold interstellar gas which
they cannot represent directly, and qualitatively reproduce the observed
$z\approx 0$ HI and H$_2$ mass functions (HIMF, H2MF), CO(1-0) luminosity
functions (COLF), and gas scaling relations versus stellar mass, specific star
formation rate, and stellar surface density $\mu_*$, with some quantitative
differences. To compare to the COLF, we apply an H$_2$-to-CO conversion factor
to the simulated galaxies based on their average molecular surface density and
metallicity, yielding substantial variations in $\alpha_{\rm CO}$ and
significant differences between models. Using this, predicted $z=0$ COLFs agree
better with data than predicted H2MFs. Out to $z\sim 2$, EAGLE's and Simba's
HIMF and COLF strongly increase, while TNG's HIMF declines and COLF evolves
slowly. EAGLE and Simba reproduce high $L_{\rm CO1-0}$ galaxies at $z\sim 1-2$
as observed, owing partly to a median $\alpha_{\rm CO}(z=2)\sim 1$ versus
$\alpha_{\rm CO}(z=0)\sim 3$. Examining \HI, H$_2$, and CO scaling relations,
their trends with $M_*$ are broadly reproduced in all models, but EAGLE yields
too little HI in green valley galaxies, TNG and Simba overproduce cold gas in
massive galaxies, and Simba overproduces molecular gas in small systems. Using
Simba variants that exclude individual AGN feedback modules, we find that
Simba's AGN jet feedback is primarily responsible by lowering cold gas contents
from $z\sim 1\to0$ by suppressing cold gas in $M_*> 10^{10}{\rm M}_\odot$
galaxies, while X-ray feedback suppresses the formation of high-$\mu_*$
systems.
