Understanding and prediction of gas-liquid interfacial dynamics is essential for several practical applications. Both sharp and diffuse interface models have been widely used for the simulation of moving interfaces under different flow and thermal conditions, providing valuable detailed information in that region. However, accurate representation of interface topological changes and phase change at interfacial regions still remains a challenge. Models for the simulation of evaporation and condensation often require the use of ad-hoc assumptions to represent mass transfer at the interface. Diffuse interface models based on the Cahn-Hilliard or Allen-Cahn equations allow to incorporate thermodynamic models, reducing the amount of ad-hoc information needed at the interface. Another approach to a diffuse representation of the interface is the Navier-Stokes Korteweg (NSK) model. Based on the theory from van der Waals, the NSK model uses the continuous varying fluid density as phase parameter and incorporates interfacial dynamics based on thermodynamic principles by modifying the stress tensor in the Navier-Stokes equations. While this formulation can deal with complex topological changes at the interface without the need for any empirical parameters, it introduces several challenges. As for any diffuse interface model, the interfacial region needs to be spatially resolved, which demands the use of very fine meshes and limits this method to small scale problems or flows near critical conditions, where the thickness of the interface becomes larger. Due to the third order derivatives of the density field appearing in the stress tensor, the system requires additional thermodynically consistent wetting boundary conditions. In addition, the NSK equations are not fully hyperbolic, which limits the application of standard compressible flow solvers. Numerical schemes such as local discontinous Galerking and relaxation models have been suggested in the literature for solving the NSK equations. This talk presents the use of the Navier-Stokes-Korteweg phase field model for the simulation of evaporation and condensation phenomena. The model is thermodynamically consistent and does not require any ad-hoc information for mass transfer at the interface separating the two phases. The least-squares spectral element method is used for the discretization of the system, using a coupled space-time formulation. The least-squares approach transforms the system into a symmetric positive definite algebraic system, favouring the use of standard solvers. A study of evaporation and condensation of droplets and thin films using this model is presented.