In this talk, we propose a Local Projection Stabilization (LPS) finite element method applied to the numerical solution of natural convection problems. Firstly, after recalling the mathematical model for which the Boussinesq approximation is employed to treat thermal coupling, a particular LPS scheme, the high-order term-by-term stabilization method, is introduced and analyzed. This method replaces the projection-stabilized structure of standard LPS methods by an interpolation-stabilized structure, which only acts on the high frequencies components of the flow. This approach gives rise to a method which may be cast in the Variational Multi-Scale (VMS) framework, and constitutes a low-cost, accurate solver for incompressible flows, despite being only weakly consistent. This method has been applied to the simulation of a high Reynolds number plane mixing layer flow, with accurate results for relatively coarse grids. Here, numerical results for the 2D problem of a buoyancy-driven airflow in a square cavity with dierentially heated side walls at high Rayleigh numbers are also given and compared with benchmark solutions. Again, a good accuracy is obtained with relatively coarse grids.