In this paper we propose and analyze, utilizing mainly tools and abstract results from Banach spaces rather than from Hilbert ones, a new fully-mixed finite element method for the stationary Boussinesq problem with temperature-dependent viscosity. More precisely, following an idea that has already been applied to the Navier-Stokes equations and to the fluid part only of our model of interest, we first incorporate the velocity gradient and the associated Bernoulli stress tensor as auxiliary unknowns. Additionally, and differently from earlier works in which either the primal or the classical dual-mixed method is employed for the heat equation, we consider here an analogue of the approach for the fluid, which consists of introducing as further variables the gradient of temperature and a vector version of the Bernoulli tensor. The resulting mixed variational formulation, which involves the aforementioned four unknowns together with the original variables given by the velocity and temperature of the fluid, is then reformulated as a fixed point equation. Next, we utilize the well-known Banach and Brouwer theorems, combined with the application of the Babuska-Brezzi theory to each independent equation, to prove, under suitable small data assumptions, the existence of a unique solution to the continuous scheme, and the existence of solution to the associated Galerkin system for a feasible choice of the corresponding finite element subspaces. Finally, we derive optimala priorierror estimates and provide several numerical results illustrating the performance of the fully-mixed scheme and confirming the theoretical rates of convergence.