In this study, we apply an artificial viscosity method to convert an unsteady level set (LS) convection equation into an unsteady LS convection-diffusion transport equation to stabilize the numerical solution of the convection term. Then a novel least-square polygonal finite element method is used to solve an unsteady LS convection-diffusion problem. The least-squares method provided good mathematical properties such as natural numerical diffusion and the positive definite symmetry of the resulting algebraic systems for the convection-diffusion and re-initialization equations. The proposed method is evaluated numerically in two different benchmark problems: a rigid body rotation of Zalesak’s disk, and a time-reversed single-vortex flow. In comparison with conventional triangular (T3) and quadrilateral (Q4) elements, polygonal elements are capable of providing greater flexibility in mesh generation for complicated problems as well as more accurate in solving the LS equations. In addition, the numerical results are also compared with the results which obtained from essentially non-oscillatory type formulations and particle LS methods. The results show that the proposed method completely matches the previously published results.