A numerical analysis was performed to study the onset of thermal instability in the thermal entrance region of horizontal rectangular ducts. Three kinds of surface thermal boundary conditions are considered. In the analysis, three dimensionless groups, namely, Rayleigh number Ra, aspect ratio-gamma, and Prandtl number Pr appear. Their corresponding influences on the local Nusselt number distributions and the onset of thermal instability are investigated. The predicted results, based on the criterion of 2% deviation of the local Nusselt number from that of Graetz theory, are in agreement with the published experimental data. Results show that the thermal instability is only slightly affected by the change of gamma under fixed Pr for the thermal condition of case 1 (the bottom wall is heated while the top wall is cooled and the side walls are adiabatic) and of case 2 (heated from below while the other walls are insulated), while for the isothermal channels (case 3), the effects of different gamma on the instability become significant. Another criterion, based on the location of the minimum local Nusselt number, is also considered in order to investigate the effects of the selection of different criteria on the prediction of the onset of thermal instability.