At first, a semi-discrete Crank–Nicolson (CN) formulation with respect to time for the non-stationary incompressible Boussinesq equations is presented. Then, a fully discrete stabilized CN mixed finite volume element (SCNMFVE) formulation based on two local Gauss integrals and parameter-free is established directly from the semi-discrete CN formulation with respect to time. Next, the error estimates for the fully discrete SCNMFVE solutions are derived by means of the standard CN mixed finite element method. Finally, some numerical experiments are presented illustrating that the numerical errors are consistent with theoretical results, the computing load for the fully discrete SCNMFVE formulation are far fewer than that for the stabilized mixed finite volume element (SMFVE) formulation with the first time accuracy, and its accumulation of truncation errors in the computational process is far lesser than that of the SMFVE formulation with the first time accuracy. Thus, the advantage of the fully discrete SCNMFVE formulation for the non-stationary incompressible Boussinesq equations is shown sufficiently.