Earth system models (ESMs) parameterize ocean surface fluxes of heat, moisture and momentum with empirical bulk flux algorithms, which introduce biases and uncertainties into simulations. We compare, for the first time, the effects of three different algorithms in both atmosphere and ocean model simulations using the Department of Energy ESM (E3SM). Flux differences between algorithms are larger in atmosphere model simulations (where wind speeds can vary) than ocean model simulations (where wind speeds are fixed by forcing data). Surface flux changes lead to global scale changes in the energy and water cycles, notably including ocean heat uptake and global mean precipitation rates. Compared to the default algorithm in E3SM, both the COARE and University of Arizona (UA) algorithms reduce global mean precipitation and top of atmosphere radiative biases. Further, the UA algorithm may slightly reduce biases in ocean meridional heat transport. We speculate that changes seen here, especially in the ocean, could be even larger in coupled simulations.