<p>This work aims to construct and analyze a virtual element method on polytopal grids to solve the fully mixed formulation of the unsteady Navier–Stokes equations coupled with the heat equation under mixed boundary conditions, commonly known as the two-dimensional nonstationary Boussinesq equations. Alongside the original thermo-fluid variables, the pseudostress and pseudoheat vector (so-called the temperature gradient) are introduced as auxiliary unknowns, driven by the growing interest in non-Newtonian flows and coupled interface problems, where stress and the thermal gradient play fundamental roles. The resulting formulation is then written equivalently as uncoupled problems thanks to a fixed-point strategy, so that the theory of differential-algebraic systems, combined with the classical Banach Theorem, is employed to establish the unique solvability of the continuous formulation. By combining the mixed virtual element approach with implicit backward Euler time integration, the problem is discretized in space-time. As a result, the proposed scheme involves nonlinear terms implicitly, and its well-posedness and stability have been established. Additionally, a convergence analysis is performed for all variables in their natural norms, demonstrating an optimal rate of convergence with respect to both the mesh size and the time step. A series of verification tests is conducted to validate the theoretical predictions and demonstrate the method’s applicability to engineering-related cases.</p>