A direct numerical simulation (DNS) of the fully developed thermal field in a two-dimensional turbulent channel flow of air was carried out. The isoflux condition was imposed on the two walls so that the local mean temperature increased linearly in the streamwise direction. With any buoyancy effect neglected, temperature was considered as a passive scalar. The computation was executed on 1,589,248 grid points by using a spectral method. The statistics obtained were root-mean-square temperature fluctuations, turbulent heat fluxes, turbulent Prandtl number, and dissipation time scales. They agreed fairly well with existing experimental and numerical simulation data. Each term in the budget equations of temperature variance, its dissipation rate, and turbulent heat fluxes was also calculated. It was found that the temperature fluctuation θ′ was closely correlated with the streamwise velocity fluctuation u′, particularly in the near-wall region. Hence, the distribution of budget terms for the streamwise and wall-normal heat fluxes, u′θ′ and v′θ′ , were very similar to those for the two Reynolds stress components, u′u′ and u′v′ .