A comparative assessment of several computational turbulence models for third-order diffusive transport terms, u iu jθ and u iθ0 2 , 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 Jrierforms equally or even better than the other more complicated ones for the scalar variance diffusion, u i 0 2. However, it is found that the computational model for the scalar flux diffusion, u i u j 0, 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.