Projection based variational multiscale (VMS) methods are a very successful technique in the numerical simulation of high Reynolds number flow problems using coarse discretizations. However, their implementation into an existing (legacy) codes can be very challenging in practice. We propose a second order variant of projection-based VMS method for non-isothermal flow problems. The method adds stabilization as a decoupled post-processing step for both velocity and temperature, and thus can be efficiently and easily used with existing codes. In this work, we propose the algorithm and give numerical results for convergence rates tests and coarse mesh simulation of Marsigli flow.