A modified thermal Bhatnagar−Gross−Krook Lattice Boltzmann (BGK-LB) model was developed to study the convection phenomenon of non-Newtonian fluids (NNFs). This model integrates the local shear rate into the equilibrium distribution function (EDF) of the flow field and keeps the relaxation time from varying with fluid viscosity by introducing an additional parameter. In addition, a modified temperature EDF was constructed for the evolution equation of the temperature field to ensure the precise recovery of the convection−diffusion equation. To validate the accuracy and effectiveness of the proposed model, numerical simulations of benchmark problems were performed. Subsequently, we investigated the natural convection of power−law (PL) fluids and examined the impact of the PL index (n = 0.7−1.3) and Rayleigh number (Ra = 103−5 × 105) on the flow and temperature fields while holding the Prandtl number (Pr = 7) constant. The obtained results indicate that, for a given value of n, the convective intensity exhibits a positive correlation with Ra, which is illustrated by the rising trend in the average Nusselt number (Nu¯) with increasing Ra. Additionally, shear-thinning fluid (n < 1) exhibited increased Nu¯ values compared to the Newtonian case, indicating an enhanced convection effect. Conversely, shear-thickening fluid (n > 1) exhibits reduced Nu¯ values, indicating weakened convective behavior.