<p>This paper proposes and studies an efficient grad-div stabilized two-grid finite element discretization method for the steady natural convection equations. In the developed method, a nonlinear natural convection problem is firstly solved on a coarse grid, and then, a Newton-linearized problem is solved on a fine grid to update the solutions. In both coarse and fine grid problems, a grad-div stabilization term is added to the variational formulation to reduce the impact of pressure on velocity error. We analyze theoretically the stability of the approximate solutions and derive the explicit dependence of errors of the approximate solutions upon Prandtl number <i>Pr</i>, Rayleigh number <i>Ra</i>, the thermal conductivity <i>κ</i>, and the stabilization parameter <i>α</i>. We finally conduct some numerical experiments to verify the theoretical analysis which show the effectiveness of the proposed method.</p>