In the context of severe nuclear accidents, the migration of corium into the reactor pressure vessel (RPV) poses significant hazards, prompting the proposal of the in-vessel melt retention (IVR) strategy, particularly the external reactor vessel cooling (ERVC) approach. Evaluating the accuracy of turbulence models within the melt pool is crucial for assessing the feasibility of IVR. However, previous studies have yet to reach a consensus about the most suitable model due to the lack of data comparison. We aim to conduct a comprehensive comparative analysis of turbulence models to evaluate their performance across a range of Rayleigh numbers, particularly under conditions relevant to IVR scenarios. Therefore, this study employs six commonly used turbulence models in computational fluid dynamics (CFD) software, ANSYS Fluent, to simulate three natural convection experiments (Kulacki-Goldstein, BALI, and LIVE-3D). The results demonstrate that the choice of turbulence model significantly impacts the accuracy of temperature and heat flux predictions within the melt pool. Although the relative temperature deviation is less than 0.1% in all the simulations of the Kulacki-Goldstein experiment, the differences among turbulence models become increasingly pronounced with rising Rayleigh numbers. Among the models tested, wall-modeled large eddy simulation (WMLES) proved the most reliable for complex geometries and high Rayleigh numbers, while the realizable k-epsilon and generalized k-omega (GEKO) models also showed consistent performance. However, the Reynolds stress model (RSM)-baseline (BSL) and detached eddy simulation (DES) models exhibited notable limitations, particularly in scenarios involving solidification and melting. These findings provide valuable guidance for selecting appropriate turbulence models in IVR-related natural convection simulations and highlight the need for further refinement to improve model accuracy across varying melt pool conditions.