Results of direct numerical simulation (DNS) for Rayleigh–Bénard convection for the Prandtl number $\hbox{\it Pr}\,{=}\,0.025$ are used to show some peculiarities of turbulent natural convection for low-Prandtl-number fluids. Simulations for this flow at sufficiently large Rayleigh numbers became feasible only recently because this flow requires the resolution of very small velocity scales and the recording of long-wave structures for the slow changes in the convective temperature field. The results are used to analyse standard turbulent heat flux models. The analysis for a model based on the Reynolds analogy indicates strong deficiencies of such turbulent heat flux models for low-Prandtl-number fluids. Turbulence models for buoyant flows which are not based on the Reynolds analogy include also the transport equation for the temperature variance $\overline{\theta^2}$. Detailed analysis of this transport equation and of the transport equation for the temperature variance dissipation rate is performed using DNS data. The results show the relevance of the turbulent diffusion terms and strong quantitative and qualitative deficiencies of standard models for turbulent diffusion of the temperature variance $\overline{\theta^2}$ and for the turbulent diffusion of the temperature variance dissipation rate $\epsilon_\theta$. Using the two-point correlation technique, statistical turbulence models for the turbulent diffusion of the temperature variance and for the turbulent diffusion of the temperature variance dissipation rate are proposed. These new models explicitly consider the molecular fluid properties. The new models reproduce the DNS results for $\hbox{\it Pr}\,{=}\,0.025$ and $\hbox{\it Pr}\,{=}\,0.71$ sufficiently well.