Direct numerical simulations (DNS) of natural convection in a vertical channel by Versteegh & Nieuwstadt (1998) are used for assessing the budget of the turbulent heat flux θui and the temperature variance θ2, and for modelling the transport equations governing these two properties. The analysis is confined to a simple fully developed situation in which the gravitational vector, as the sole driving force, is perpendicular to the only non-zero component of the mean temperature gradient. Despite its simplicity, the flow displays many interesting features and represents a generic case of the interaction of buoyancy-driven turbulent temperature and velocity fields. The paper discusses the near-wall variation of the second moments and their budgets, as well as possible scaling of θui and θ 2 both in the near-wall region and away from the wall. Various proposals for the Reynolds-averaged modelling are analysed and new models are proposed for these two transport equations using the term-by- term approach. An a priori test (using the DNS data for properties other than θui and θ 2) reproduced very well all terms in the transport equations, as well as their near-wall behaviours and wall limits, without the use of any wall-topology-dependent parameters. The computational effort is still comparable to that for the ‘basic model’. The new term-by-term model of the θui and θ 2 equations was then used for a full simulation in conjunction with a low-Reynolds-number second-moment velocity closure, which was earlier found to reproduce satisfactorily a variety of isothermal wall flows. Despite excellent term-by-term reproduction of thermal turbulence, the predictions with the full model show less satisfactory agreement with the DNS data than a priori validation, indicating a further need for improvement of the modelling of buoyancy effects on mechanical turbulence.