The numerical simulation of turbulent heat transfer at high molecular Prandtl numbers is still strongly challenged by disparate thicknesses of the viscous sublayer for velocity and enthalpy. The model assumption for the prescription of the thermal boundary conditions in Reynolds Averaged Navier-Stokes (RANS) type simulations becomes increasingly questionable for high molecular Prandtl numbers, where the Reynolds analogy between the momentum and the heat transfer breaks down.The present work addresses this issue by performing Direct Numerical Simulation (DNS) for relatively high molecular Prandtl number, up to Pr=50, including temperature-dependent fluid properties. The variation of all flow quantities, due to varying fluid properties is extensively examined by a comparison against results from DNS with constant fluid properties.All DNS results are further validated against measurements, obtained on a test facility, specially designed and constructed for this purpose. The agreement of the simulated skin-friction coefficient and Nusselt numbers with the experiments is excellent, indicating the high reliability of the presented DNS results.This DNS database was utilized for the evaluation and further development of two popular modeling approaches for RANS-type thermal boundary condition, the so-called P-function and a simple two-layer model. Based on the DNS results, various modifications for either approach are proposed, providing in particular appropriate submodels for the turbulent Prandtl number, the eddy viscosity, and incorporating an analytical solution for the ratio of the total shear stress and the total heat flux. The resulting modified formulations are capable of reproducing very accurately the DNS results as well as the experimental data which covers also a Reynolds numbers not amenable to DNS.