  We investigate the evolution of the ortho-to-para ratio of overall (gas +
ice) H$_2$ via the nuclear spin conversion on grain surfaces coated with water
ice under physical conditions that are relevant to star- and planet-forming
regions. We utilize the rate equation model that considers adsorption of
gaseous H$_2$ on grain surfaces which have a variety of binding sites with a
different potential energy depth, thermal hopping, desorption, and the nuclear
spin conversion of adsorbed H$_2$. It is found that the spin conversion
efficiency depends on the H$_2$ gas density and the surface temperature. As a
general trend, enhanced H$_2$ gas density reduces the efficiency, while the
temperature dependence is not monotonic; there is a critical surface
temperature at which the efficiency is the maximum. At low temperatures, the
exchange of gaseous and icy H$_2$ is inefficient (i.e., adsorbed H$_2$ does not
desorb and hinders another gaseous H$_2$ to be adsorbed), while at warm
temperatures, the residence time of H$_2$ on surfaces is too short for the spin
conversion. Additionally, the spin conversion becomes more efficient with
lowering the activation barriers for thermal hopping. We discuss whether the
spin conversion on surfaces can dominate over that in the gas-phase in star-
and planet-forming regions. Finally, we establish a simple but accurate way to
implement the H$_2$ spin conversion on grain surfaces in existing astrochemical
models.
