  This paper considers a modular grad-div stabilization method for
approximating solutions of the time-dependent Boussinesq model of
non-isothermal flows. The proposed method adds a minimally intrusive step to an
existing Boussinesq code, with the key idea being that the penalization of the
divergence errors, is only in the extra step (i.e. nothing is added to the
original equations). The paper provides a full mathematical analysis by proving
unconditional stability and optimal convergence of the methods considered.
Numerical experiments confirm theoretical findings, and show that the
algorithms have a similar positive effect as the usual grad-div stabilization.
