In this paper we complement the study of a new mixed finite element scheme, allowing conservation of momentum and thermal energy, for the Boussinesq model describing natural convection and derive a reliable and efficient residual-based a posteriori error estimator for the corresponding Galerkin scheme in two and three dimensions. More precisely, by extending standard techniques commonly used on Hilbert spaces to the case of Banach spaces, such as local estimates, suitable Helmholtz decompositions and the local approximation properties of the Clément and Raviart–Thomas operators, we derive the aforementioned a posteriori error estimator on arbitrary (convex or non-convex) polygonal and polyhedral regions. In turn, inverse inequalities, the localization technique based on bubble functions, and known results from previous works, are employed to prove the local efficiency of the proposed a posteriori error estimator. Finally, to illustrate the performance of the adaptive algorithm based on the proposed a posteriori error indicator and to corroborate the theoretical results, we provide some numerical examples.