In this article, an operator‐splitting (or a fractional‐step) finite element method is proposed for the numerical solution of the time‐dependent natural convection problem. We first analyze the time‐discrete fractional‐step scheme, which decouples the considered system into three subproblems, and each subproblem can be solved more easily than the original one. Under some mild regularity assumptions, the stability and error estimates for the velocity and the temperature are established rigorously. In addition, a full discrete fractional‐step scheme based on the time‐discrete scheme is proposed and the rigorous error analysis is presented. We derive the temporal–spatial error estimates of scriptO(normalΔt+h)$$ \mathcal{O}\left(\Delta t+h\right) $$ for the velocity and the temperature in the discrete space l2()H1∩l∞()L2$$ {l}^2\left({H}^1\right)\cap {l}^{\infty}\left({L}^2\right) $$ under the constraint Δt ≥ Ch. Numerical experiments are performed to confirm the theoretical predictions and demonstrate the efficiency of the method.