We present a novel numerical method to solve the incompressible Navier-Stokes equations for two-phase flows with phase change, using a one-fluid approach. Separate phases are tracked using a geometric Volume-Of-Fluid (VOF) method with piecewise linear interface construction (PLIC). Thermal energy advection is treated in conservative form and the geometric calculation of VOF fluxes at computational cell boundaries is used consistently to calculate the fluxes of heat capacity. The phase boundary is treated as sharp (infinitely thin), which leads to a discontinuity in the velocity field across the interface in the presence of phase change. The numerical difficulty of this jump is accommodated with the introduction of a novel two-step VOF advection scheme. The method has been implemented in the open source code PARIS and is validated using well-known test cases. These include an evaporating circular droplet in microgravity (2D), the Stefan problem and a 3D bubble in superheated liquid. The accuracy shown in the results were encouraging. The 2D evaporating droplet showed excellent prediction of the droplet volume evolution as well as preservation of its circular shape. A relative error of less than 1% was achieved for the Stefan problem case, using water properties at atmospheric conditions. For the final radius of the bubble in superheated liquid at a Jacob number of 0.5, a relative error of less than 6% was obtained on the coarsest grid, with less than 1% on the finest.