The objective of this work is to describe a variational multiscale finite element approximation for the incompressible Navier-Stokes equations using the Boussinesq approximation to model thermal coupling. The main feature of the formulation in contrast to other stabilized methods is that we consider the subscales 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, we keep the effect of the subscales both in the nonlinear convective terms of the momentum and temperature equations and, if required, in the thermal coupling term of the momentum equation. This strategy allows us to approach the problem of dealing with thermal turbulence from a strictly numerical point of view and discuss important issues, such as the relationship between the turbulent mechanical dissipation and the turbulent thermal dissipation.