Purpose - The purpose of this paper is to describe a variational multiscale finite element approximation for the incompressible Navier-Stokes equations using the Boussinesq approximation to model thermal coupling. Design/methodology/approach - The main feature of the formulation, in contrast to other stabilized methods, is that the subscales are considered as transient and orthogonal to the finite element space. These subscales are solution of a differential equation in time that needs to be integrated. Likewise, the effect of the subscales is kept, both in the nonlinear convective terms of the momentum and temperature equations and, if required, in the thermal coupling term of the momentum equation. Findings - This strategy allows the approaching of the problem of dealing with thermal turbulence from a strictly numerical point of view and discussion important issues, such as the relationship between the turbulent mechanical dissipation and the turbulent thermal dissipation. Originality/value - The treatment of thermal turbulence from a strictly numerical point of view is the main originality of the work.