In this paper, we report on the efficacy of four different thermal boundary models for Wall-Modelled Large Eddy Simulations (WMLES) of turbulent natural convection. Our test cases consist of Rayleigh-Bénard convection of liquid water at two Rayleigh numbers, \(Ra =1.35{\times }10^8\) and \(Ra =10^9\) , respectively. Two configurations are examined, namely, convection in a box and in a cavity; the latter one involving a free-slip top boundary. For these test cases, the numerical results obtained via WMLES with the thermal boundary models are compared with those of Wall-Resolved Large-Eddy Simulations. According to our comparative studies, a particular version of the so-called Kays & Crawford model provides the most accurate predictions, at least for the test cases considered herein. Additionally, in this paper, we report on WMLES of turbulent convection at a higher Rayleigh number, \(Ra =5{\times }10^9\) , with the aforementioned model. For this case, we analyse herein the flow structure and present results for first and second-order statistics of the flow.