  The purpose of this article is to study the convergence of a low order finite
element approximation for a natural convection problem. We prove that the
discretization based on P1 polynomials for every variable (velocity, pressure
and temperature) is well-posed if used with a penalty term in the divergence
equation, to compensate the loss of an inf-sup condition. With mild assumptions
on the pressure regularity, we recover convergence for the
Navier-Stokes-Boussinesq system, provided the penalty term is chosen in
accordance with the mesh size. We express conditions to obtain optimal order of
convergence. We illustrate theoretical convergence results with extensive
examples. The computational cost that can be saved by this approach is also
assessed.
