The primary objective of the present study is to investigate the influence of magnetohydrodynamic (MHD) flow and heat transfer behavior of Jeffrey fluid under natural convection within a square cavity filled with a permeable matrix. This investigation is significant because enhancing heat transfer capabilities in systems such as nuclear reactor cooling is crucial for ensuring efficient thermal management. The cavity is configured with cold vertical walls, an adiabatic top surface, and a heated bottom surface, while a constant vertical magnetic field is applied at the left wall. The momentum transfer in the permeable matrix is modelled using the Darcy–Forchheimer approach, and the Galerkin finite element method (GFEM) is employed within COMSOL Multiphysics 6.1 to solve the governing equations. The study examines a range of Rayleigh numbers (10³ ≤ Ra ≤ 106), Darcy numbers (10-5 ≤ Da ≤ 10-3), and Hartmann numbers (10 ≤ Ha ≤ 40), providing a detailed analysis of the Nusselt number, velocity distribution, isocontours, isotherms, temperature profiles, and stream functions. Key findings of the study reveal that as the Hartmann number increases, the velocity distribution exhibits a monotonic rise which indicating the strong influence of the magnetic field on flow dynamics. Numerical results of the study demonstrate that with an increase in the Hartmann number (Ha) from 10 to 40, the average Nusselt number on the hot wall decreases from 13.645 to 12.380 at a Rayleigh number (Ra) of 106 indicate a reduction in heat transfer efficiency due to the damping effect of the magnetic field. For lower Rayleigh numbers (Ra = 103) the Nusselt number remains nearly constant around 5.728 across varying Hartmann numbers which shows that the magnetic field’s impact is less significant under weaker convective conditions. The results of the study show a high degree of consistency with previous studies, demonstrating the robustness of the numerical approach. This work advances the understanding of MHD natural convection with Jeffrey fluids by offering specific, quantitative insights that go beyond previous literature, particularly in the context of optimizing heat transfer in engineering applications. The novelty of present findings are particularly relevant to geophysical applications, such as modeling the movement of magma in volcanic cavities, as well as industrial processes like polymer mixing.