The present study numerically investigates the performance of a range of RANS turbulence closures with different near-wall treatments in the prediction of natural convection flows. Comparisons are performed of both 2-D and 3-D steady and time-dependent computations in differentially heated cavities, which involve very simple geometries, but give rise to complex flow physics of relevance to cooling applications, including safety features of nuclear reactors. The turbulence closures assessed include eddy-viscosity and second-moment closures and involve low-Reynolds-number versions (LRN) which fully resolve the viscous sublayer, as well as high-Reynolds-number versions (HRN) which adopt wall function treatments to prescribe boundary conditions for the wall shear stress and heat flux. In the 2-D computations of the tall rectangular cavity, there are only small deviations in the predictions of different models, and these are mainly confined to the local Nusselt number comparisons. Among the LRN models, the second-moment EBRSM is in good agreement with the data, while the LRN eddy-viscosity models under-estimate the Nusselt number levels. The HRN models, when used with the standard log-law-based wall function, over-estimate the Nusselt umber. During the course of the present investigation, a new, more general, version of the Analytical Wall Function (AWF) has been developed and introduced to the open-source CFD code utilised. Here, this more general numerical wall treatment brings the HRN k − ε predictions of the Nusselt number to very close agreement with the experimental data. In the 2-D square cavity computations, the simulations focus on a high Rayleigh number (Ra) of 1011, for which DNS data has recently become available. At such high Ra value all the models are challenged due to the unsteadiness of the flow and the presence of a largely stationary and laminar core, and thin turbulent boundary layers. Assessment of different 3-D configurations of the same square cavity demonstrates that the effect of the turbulent mixing in the third (z) direction is significant, but can be taken into account by introducing a cavity depth of D = 0.15H, matching that used in the DNS study. The deviations in the predictions of the different models are stronger than those observed in the corresponding rectangular cavity comparisons, especially for the local Nusselt number distribution. As is also the case in the rectangular cavity, the eddy-viscosity models under-estimate the Nusselt number, while the EBRSM with the more elaborate turbulent heat flux model once more is closest to the DNS data. All HRN models with the log-law wall function severely underestimate the Nusselt number levels. Introduction of the buoyancy-adapted AWF leads to substantial predictive improvements. The HRN AWF approach is thus shown to be a reliable and cost-effective modelling strategy for natural convection flows.