This study presents a modular sparse grad-div stabilization method for solving the Boussinesq equations. Unlike the usual grad-div stabilization which produces fully coupled block matrices, the proposed stabilization method produces block upper triangular matrices. Thus, the proposed method is more attractive in terms of both its computational cost and solution accuracy. We provide unconditional stability results for velocity and temperature. Two numerical experiments are performed to demonstrate the efficiency and accuracy of the method.