The objectives of this study are to reconstruct a turbulence model of both the large Eddy simulation (LES) and the Reynolds-averaged Navier–Stokes simulation (RANS) which can predict wind synopsis in various thermally stratified turbulent boundary layers over any obstacles. Hence, the direct numerical simulation (DNS) of various thermally stratified turbulent boundary layers with/without forward-step, two-dimensional block, or two-dimensional hill is carried out in order to obtain detailed turbulent statistics for the construction of a database for the evaluation of a turbulence model. Also, DNS clearly reveals the characteristics of various thermally stratified turbulent boundary layers with/without forward-step, two-dimensional block, or two-dimensional hill. The turbulence models employed in LES and RANS are evaluated using the DNS database we obtained. In the LES, an evaluated turbulence model gives proper predictions, but the quantitative agreement of Reynolds shear stress with DNS results is difficult to predict. On the other hand, the nonlinear eddy diffusivity turbulence models for Reynolds stress and turbulent heat flux are also evaluated using DNS results of various thermally stratified turbulent boundary layers over a forward-step in which the turbulence models are evaluated using an a priori method. Although the evaluated models do not make it easy to properly predict the Reynolds shear stresses in all cases, the turbulent heat fluxes can be qualitatively predicted by the nonlinear eddy diffusivity for a heat turbulence model. Therefore, the turbulence models of LES and RANS should be improved in order to adequately predict various thermally stratified turbulent boundary layers over an obstacle.