The effect of Prandtl number (Pr) on the multiplicity of the steady states (flow and temperature maps) is investigated for the natural convection of Newtonian fluids within a square enclosure involving hot bottom wall, cold side walls and adiabatic top wall. The analysis is based on the steady state governing equations, which are solved by penalty Galerkin finite element method and Newton–Raphson solution scheme. A novel perturbation technique has been developed, which is applied in conjunction with an in-house continuation scheme for the initiation of various symmetric as well as asymmetric steady states. The current computation strategy involving perturbation and continuation schemes leads to a rich bifurcation diagram in the parameter space of Pr (0.01⩽Pr⩽10) for Ra (Rayleigh number) =106. It is found that the occurrence of multiple steady states depends on Pr. Only one symmetric solution is obtained in the high–Pr regime (Pr>0.42), while a highly populated bifurcation diagram consisting of 6 symmetric branches, 18 asymmetric branches and 4 isolas is obtained in the low–Pr regime (especially for Pr≲0.05). The solution branches correspond to a wide spectrum of flow structures involving single cell, multiple cells, reverse-cells, corner-cells, horizontally stacked cells, vertically stacked cells etc., which are reported for the first time. These flow structures result in a wide range of Nusselt numbers (average heat transfer rates) at a given Pr, especially for the low–Pr fluids (such as liquid metals). The flow structures are critically analyzed in terms of the location of the flow intensification and its impact on the overall heat transfer rates. It has been shown that the overall heat transfer rates for the low–Pr fluids can be significantly enhanced (80 –100% in few cases) based on corner-cells or single cell driven flow. © 2021 Elsevier Ltd