In this work we present and analyse a new fully mixed finite element method for the nonlinear problem given by the coupling of the Darcy and heat equations. Besides the velocity, pressure and temperature variables of the fluid, our approach is based on the introduction of the pseudoheat flux as a further unknown. As a consequence of it, and due to the convective term involving the velocity and the temperature, we arrive at saddle point-type schemes in Banach spaces for both equations. In particular, and as suggested by the solvability of a related Neumann problem to be employed in the analysis, we need to make convenient choices of the Lebesgue and H(div)-type spaces to which the unknowns and test functions belong. The resulting coupled formulation is then written equivalently as a fixed-point operator, so that the classical Banach theorem, combined with the corresponding Babuska-Brezzi theory, the Banach-Necas-Babuska theorem, suitable operators mapping Lebesgue spaces into themselves, regularity assumptions and the aforementioned Neumann problem, are employed to establish the unique solvability of the continuous formulation. Under standard hypotheses satisfied by generic finite element subspaces, the associated Galerkin scheme is analysed similarly and the Brouwer theorem yields existence of a solution. The respective a priori error analysis is also derived. Then, Raviart-Thomas elements of order k >= 0 for the pseudoheat and the velocity and discontinuous piecewise polynomials of degree <= k for the pressure and the temperature are shown to satisfy those hypotheses in the two-dimensional case. Several numerical examples illustrating the performance and convergence of the method are reported, including an application into the equivalent problem of miscible displacement in porous media.