The classical thermodynamic properties of water clusters modeled by the TIP4P potential with and without the inclusion of molecular hydrogen are calculated using Monte Carlo methods. The temperature-dependent heat capacity curves of the dimer and the smaller pure water clusters having ring structures show low temperature anomalies arising from the onset of a transition from librational motion to free rotational motion. Pure water clusters having cage structures display heat capacity anomalies characteristic of "melting" phase changes. The addition of molecular hydrogen to the water clusters has little impact on the structure of the water core, but low temperature heat capacity anomalies are observed that arise from the onset of hindered to free rotational and translational motion of the hydrogen molecule with respect to the skeletal moieties of the core. The Gibbs free energy changes associated with the growth of pure water clusters and the addition of molecular hydrogen to the water clusters are determined using a combination of stare-integration methods and an application of the Gibbs-Helmholtz equation. For the temperature integration of the Gibbs-Helmholtz equation, a quadrature method is introduced that avoids numerical difficulties arising from singularities in the integrand at low temperatures. For the growth of pure water clusters, the fine structure of the enthaply and the low-temperature Gibbs free energy as a function of cluster size is rationalized using an ansatz of molecular dipole orientations. For the addition of molecular hydrogen to water clusters, the Gibbs free energy change is a virtually flat function of cluster size showing no fine structure.