This paper proposes a class of semi-discrete and fully discrete weak Galerkin finite element methods for unsteady natural convection problems in two and three dimensions. In the space discretization, the methods use piecewise polynomials of degrees $k,$ $k-1,$ and $k$ $(k\geq 1)$ for the velocity, pressure and temperature approximations in the interior of elements, respectively, and piecewise polynomials of degree $k$ for the numerical traces of velocity, pressure and temperature on the interfaces of elements. In the temporal discretization of the fully discrete method, the backward Euler difference scheme is adopted. The semi-discrete and fully discrete methods yield globally divergence-free velocity solutions. Wellposedness of the semi-discrete scheme is established and a priori error estimates are derived for both the semi-discrete and fully discrete schemes. Numerical experiments demonstrate the robustness and efficiency of the methods.