In geothermal modelling, a stable method for computing the thermodynamic properties of water – density, temperature and viscosity – is essential. Instability can lead to convergence issues within the flow simulator. For water, the most common way of computing those properties is through using either Gibbs or Helmholtz equations depending on the pressure (𝑝) and enthalpy (ℎ) region the fluid is in. However, the boundaries between the different regions present some discontinuities in the values of fluid properties, which can cause slow convergence of the flow simulator. In the past, these discontinuities have been addressed through the introduction of blending zones between regions to ensure a smooth transition. Our work aims to implement an efficient method of computing fluid properties that builds in property and derivative continuity across the entire 𝑝-ℎ domain. We do this by using bivariate splines with optimised, non-uniform knot locations. In addition to increased stability of the computation, this grants easy access to the derivatives of fluid properties with respect to pressure and enthalpy. This is important for some methods of reservoir simulation and optimization. The methodology presented in this paper is as follows. First, we assemble data related to the water properties we are interested in. Second, we introduce an algorithm capable of building a bivariate spline approximation that captures the complexity of our datasets while remaining sufficiently smooth. Third, we perform a quality check on the relative errors, partial derivatives and continuity of our surfaces. Finally, we implement our thermodynamic formulation into a flow simulator in order to assess its performance and precision.