In-vessel thermal stratification analysis was carried out using a multidimensional thermohydraulic analysis code, in which a higher-order finite difference scheme was applied to the convection terms. Discussions centred on the buoyancy modelling in the vicinity of the stratification interface through comparisons between experiment and calculation. Computational results were obtained from the following three turbulence models: (i) the k−ϵ model with a constant turbulent Prandtl number Pr t, (ii) the k−ϵ model with the turbulent Prandtl number being dependent on the local Richardson number Ri, and (iii) the algebraic stress model. Numerical analysis of the stratification phenomena using the higher-order scheme showed that, in general, the modelling of the buoyancy terms appearing in the turbulence transport equations was the most important key to successful results. When the k−ϵ model was used, it was pointed out that a dependence on the local Richardson number must be carefully included in the turbulent Prandtl number. In this case, however, the range of applicability was limited to the phenomena observed in the water system in general because the model was constructed and calibrated for water experiments. Overall it was found that the calculated stratification interface rise agreed well with experimental results in water and sodium insofar as the algebraic stress model was utilized. As a conclusion, in predicting the behaviour of the thermal stratification phenomena in liquid metal cooled reactors, the coupled use of the higher-order difference scheme and the algebraic stress model was most appropriate and recommended.