We present a multi-physics model for the approximation of the coupled system formed by the temperature-dependent Navier-Stokes equations with free surfaces. The main application is the industrial process of shallow laser surface melting (SLSM), for laser polishing of metal surfaces. We consider incompressible flow equations with solidification, and we model the laser source through physically-consistent boundary conditions. We incorporate Marangoni effects in the surface tension model to drive internal motion in the liquid metal. The numerical method relies on an operator splitting strategy and a two-grid approach. A proof of concept of the numerical model is achieved through a static laser melting process.