A comparative assessment of several computational turbulence models for third-order diffusive transport terms, u(i)u(j)thetaBAR and u(i)theta2BAR, in second-order closure equations has been carried out by applying the models to various non-isothermal turbulent flows. The second-order quantities appearing in the models are adopted from directly measured values. The models tested in the present study are: conventional simple gradient model, eddy-damped quasi-normal approximation model, and Weinstock's theoretical model which is derived by formally integrating the Navier-Stokes equation [J. Weinstock, J. Fluid Mech. 202, 319-338 (1989)]. It is rather a surprise to find that the simple gradient model performs equally or even better than the other more complicated ones for the scalar variance diffusion, u(i)theta2BAR. However, it is found that the computational model for the scalar flux diffusion, u(i)u(j)thetaBAR, must include the shear-gradient contribution in addition to the simple gradient model. Moreover, a buoyancy correction method is proposed to take into account the buoyancy effect in the gradient-type models.