The residual Gibbs energy surface, i.e. the residual Gibbs energy G res (P, T ) as a function of temperature and pressure, contains all the information with respect to changes in temperature and pressure. Specifically, all other property surfaces as a function of T and P are computable by derivatives of temperature and pressure of the residual Gibbs energy. So one might intuitively think that parameterizing with Gibbs energy alone should be capable of finding the most accurate water model with respect to all thermodynamic properties. However, our experience shows (as detailed in Section 1.2 of the SI), we need a very finely spaced grid in T, P space to accurately reproduce the experimental Gibbs energy surface with the all the gradients and curvatures properly captured. This is because the property surfaces can bend, twist and warp differently at different points on the T, P surface in many different ways as the parameters vary. A grid sufficiently closely spaced to make all derivatives smooth makes optimization process computationally prohibitive. To create smooth surfaces, we thus constrain the magnitude of derivatives of the surface, through the bounds on the errors in the physical properties that depend on the derivatives. Specifically, the variations in the Gibbs energy surface are related to the the gradients and curvatures of the Gibbs energy surface which in turn are related to the first, second and the mixed derivative thermodynamic properties such as density, enthalpy, heat capacity and others. One of our main findings is that adding constraints on these properties can therefore significantly improve the regularity of the parameterization.