A technique, based on the Lund et al. (1998) [1] approach, is introduced in order to numerically prescribe time-dependent turbulent inflow conditions in spatially-developing thermal boundary layers. A major difference with Lund's approach is that the new multi-scale approach considers different scaling laws for the inner and outer parts of the boundary layer. Direct numerical simulations (DNS) are performed for incompressible zero (ZPG) and adverse (APG) pressure gradient flows. To the best of our knowledge, the present DNS in ZPG flows at a momentum thickness Reynolds numbers (Reo) of 2300 represents the evolving thermal boundary layer simulations at the highest Reo in the turbulence community. The temperature is treated as a passive scalar with isothermal walls as a boundary condition and a molecular Prandtl number of 0.71. The predicted Stanton number shows fairly good agreement with empirical correlations, and experimental and numerical data from the literature. Moreover, the influences of the Reynolds number and the APG strength on thermal parameters are also examined. Furthermore, the budget of the temperature variance shows a significant increase of production, turbulent diffusion, and dissipation in the buffer layer at higher Reynolds numbers. The main effects of adverse pressure gradients on the temperature field are manifested by a decreasing trend of thermal fluctuations but an increasing trend of the wall-normal turbulent heat fluxes when normalized by wall units as the APG strength is augmented. (C) 2012 Elsevier Ltd. All rights reserved.