This paper presents stabilized finite element simulations of unsteady natural and mixed heat convection phenomena occurring in square enclosures. It is assumed that the enclosures are subject to intense uniform magnetic fields applied perpendicular to the vertical walls. It is well known that the classical Galerkin finite element method (GFEM) is prone to yielding nonphysical oscillations when simulating convection-dominated flows. Therefore, in order to handle such flows occurring at high Hartmann (0≤Ha≤100\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$0 \\le {\\rm Ha} \\le 100$$\\end{document}) and Rayleigh (103≤Ra≤108\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$10^{3} \\le {\\rm Ra} \\le 10^{8}$$\\end{document}) numbers, the stabilization of the GFEM is accomplished with the streamline-upwind/Petrov–Galerkin and pressure-stabilizing/Petrov–Galerkin formulations. Temporal discretizations are performed with the backward Euler formula. At each time step, the nonlinear systems are linearized with the Newton–Raphson (N–R) method, and the linearized systems are solved directly using the sparse lower–upper decomposition. The incompressible flow solvers are developed in-house, run in parallel within the FEniCS ecosystem, and tested on a comprehensive set of numerical experiments with various thermal boundary conditions. Numerical simulations and comparisons reveal that the proposed formulation performs quite well at high Hartmann and Rayleigh numbers without exhibiting any significant localized or globally spread numerical instabilities. Moreover, it achieves this by using uniform meshes and only linear elements, making the formulation much more computationally efficient. Mesh convergence analyses are also performed to demonstrate the reliability of the numerical results obtained.