A thermal lattice Boltzmann method based on the BGK model has been used to simulate high Rayleigh number natural convection in a square cavity. The model uses the double populations approach to simulate hydrodynamic and thermal fields. The traditional lattice Boltzmann method on a uniform grid has unreasonably high grid requirements at higher Rayleigh numbers which renders the method impractical. In this work, the interpolation supplemented lattice Boltzmann method has been utilized. This is shown to be effective even at high Rayleigh numbers. Numerical results are presented for natural convection in a square cavity with insulated horizontal walls and isothermal vertical walls maintained at different temperatures. Very fine grids (wall y+ < 0.3) have been used for the higher Rayleigh number simulations. A universal structure is shown to exist in the mean velocity turbulent boundary layer profile for y+ < 10. This agrees extremely well with previously reported experimental data. The numerical results (for Rayleigh numbers up to 1010) are in very good agreement with the benchmark results available in the literature. The highlight of the calculations is that no turbulence model has been employed. © 2005 Elsevier Ltd. All rights reserved.