  This paper proposes and analyzes a class of weak Galerkin (WG) finite element
methods for stationary natural convection problems in two and three dimensions.
We use piecewise polynomials of degrees k, k-1, and k(k>=1) for the velocity,
pressure, and temperature approximations in the interior of elements,
respectively, and piecewise polynomials of degrees l, k, l(l = k-1,k) for the
numerical traces of velocity, pressure and temperature on the interfaces of
elements. The methods yield globally divergence-free velocity solutions.
Well-posedness of the discrete scheme is established, optimal a priori error
estimates are derived, and an unconditionally convergent iteration algorithm is
presented. Numerical experiments confirm the theoretical results and show the
robustness of the methods with respect to Rayleigh number.
