Natural convection heat transfer phenomena in nanofluid-filled square cavities with various temperature boundary conditions are studied computationally. From electronic cooling to building ventilation systems, such phenomena have numerous practical applications, and accurate simulations are crucial for developing new designs. Towards that end, the Navier–Stokes equations of incompressible flows are considered with thermal coupling. The base fluid is pure water, the nanoparticles are copper (Cu), cupric oxide (CuO), or aluminum oxide (Al2O3), and the nanofluids are assumed to be homogeneous. It is well known that, in the standard finite element method framework, inappropriate selection of interpolation functions, e.g., (bi-)linear equal-order-interpolation velocity-pressure (e.g., and ) elements, yields nonphysical oscillations in the flow field for simulating incompressible flows, particularly for high Rayleigh numbers. In this study, in order to overcome such numerical instability issues, the streamline-upwind/Petrov–Galerkin (SUPG) and pressure-stabilizing/Petrov–Galerkin (PSPG) finite element formulations are utilized. The SUPG/PSPG-stabilized (SUPS) formulation is also enhanced with the least-squares on incompressibility constraint (LSIC). A comprehensive set of numerical test computations is considered for the values of the Rayleigh numbers ranging from to and a broad range of volume fractions of nanoparticles from to . Incompressible flow solvers are developed in-house and executed in parallel. Numerical simulations and comparisons with reported studies reveal that the proposed formulation performs quite well even at high Rayleigh numbers, and it does not exhibit any significant local or globally spread numerical instabilities. It is also noted that this is achieved without employing any adaptive mesh strategies and using only linear and equal-order interpolation functions, which in turn saves computational time.