In this paper we focus on the analysis of a mixed finite element method for a class of natural convection problems in two dimensions. More precisely, we consider a system based on the coupling of the steady-state equations of momentum (Navier-Stokes) and thermal energy by means of the Boussinesq approximation (coined the Boussinesq problem), where we also take into account a temperature dependence of the viscosity of the fluid. The construction of this finite element method begins with the introduction of the pseudostress and vorticity tensors, and a mixed formulation for the momentum equations, which is augmented with Galerkin-type terms, in order to deal with the non-linearity of these equations and the convective term in the energy equation, where a primal formulation is considered. The prescribed temperature on the boundary becomes an essential condition, which is weakly imposed, leading us to the definition of the normal heat flux through the boundary as a Lagrange multiplier. We show that this highly coupled problem can be uncoupled and analysed as a fixed-point problem, where Banach and Brouwer theorems will help us to provide sufficient conditions to ensure well-posedness of the problems arising from the continuous and discrete formulations, along with several applications of continuous injections guaranteed by the Rellich-Kondrachov theorem. Finally, we show some numerical results to illustrate the performance of this finite element method, as well as to prove the associated rates of convergence.