A numerical solution for the turbulent convective heat transfer affected by buoyancy force is presented. Recent numerical investigations of the problem have been carried out by the two-equation model derived from the truncated Reynolds stress closure which includes buoyancy terms. In the past, the turbulent time scale ratio R of temperature relative to velocity is assumed to be constant, which is not appropriate when the buoyancy effect is dominant. For the present analysis, R is considered as variable and a four-equation model is set up adding two equations of the mean square temperature variance $\overline{Θ}^2$ and its rate of dissipation $ε_Θ$, and is applied for the turbulent boundary layer over a heated or cooled horizontal flat plate. It is found that, R varies widely across the boundary layer and four-equation model yields better results for the mean temperature profiles and the surface heat flux for flows with large buoyancy effect.