ABSTRACT We present a comparison of galaxy atomic and molecular gas properties in three recent cosmological hydrodynamic simulations, namely SIMBA, EAGLE, and IllustrisTNG, versus observations from z ∼ 0 to 2. These simulations all rely on similar subresolution prescriptions to model cold interstellar gas that they cannot represent directly, and qualitatively reproduce the observed z ≈ 0 H i and H2 mass functions (HIMFs and H2MFs, respectively), CO(1–0) luminosity functions (COLFs), and gas scaling relations versus stellar mass, specific star formation rate, and stellar surface density μ*, with some quantitative differences. To compare to the COLF, we apply an H2-to-CO conversion factor to the simulated galaxies based on their average molecular surface density and metallicity, yielding substantial variations in αCO and significant differences between models. Using this, predicted z = 0 COLFs agree better with data than predicted H2MFs. Out to z ∼ 2, EAGLE’s and SIMBA’s HIMFs and COLFs strongly increase, while IllustrisTNG’s HIMF declines and COLF evolves slowly. EAGLE and simba reproduce high-LCO(1–0) galaxies at z ∼ 1–2 as observed, owing partly to a median αCO(z = 2) ∼ 1 versus αCO(z = 0) ∼ 3. Examining H i, H2, and CO scaling relations, their trends with M* are broadly reproduced in all models, but EAGLE yields too little H i in green valley galaxies, IllustrisTNG and SIMBA overproduce cold gas in massive galaxies, and SIMBA overproduces molecular gas in small systems. Using SIMBA variants that exclude individual active galactic nucleus (AGN) feedback modules, we find that SIMBA’s AGN jet feedback is primarily responsible by lowering cold gas contents from z ∼ 1 → 0 by suppressing cold gas in $M_*\gtrsim 10^{10}{\rm \,M}_⊙$ galaxies, while X-ray feedback suppresses the formation of high-μ* systems.