In this paper, we propose and study a numerical scheme for the Navier-Stokes-Fourier system derived and studied by the authors in [12]. This system models hydrostatic free surface flows with density variations depending on temperature or salinity. We show that the finite volume/finite element scheme presented – based on a layer averaged formulation of the model – is well-balanced with regards to the steady state of the lake at rest and preserves the nonnegativity of the water height. A maximum principle on the density is also proved as well as a discrete entropy inequality (when the thermal and viscous effects are neglected). Some numerical validations are finally shown with comparisons to 3D analytical solutions and experiments.