In this paper, we propose 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. They 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. Apart from presenting the main properties of the formulation, we also discuss some computational aspects such as the linearization strategy or the way to integrate in time the equation for the subscales.