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 <(theta u(i))over bar> and the temperature variance <(theta(2))over bar>, 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 <(theta u(i))over bar> and <(theta(2))over bar> 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 <(theta u(i))over bar> and <(theta(2))over bar>) 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 <(theta u(i))over bar> and <(theta(2))over bar> 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.