In this paper, we adapt a modified internal energy non-equilibrium first-order extrapolation thermal boundary condition to the thermal lattice Boltzmann model (TLBM). This model is the double populations approach to simulate hydrodynamic and thermal fields. The bounce-back boundary condition which is a traditional boundary condition of lattice Boltzmann method has only a first order in numerical accuracy at the boundary and numerical instability. A non-equilibrium first-order extrapolation boundary condition has been verified to be of better numerical stability than the bounce-back boundary condition and this boundary condition is proved to be of second-order accuracy for the flat boundaries. The two-dimensional natural convection flow in a square cavity with Pr=0.71 and various Rayleigh numbers are simulated. The results are found to be in good agreement with those of previous studies.