Abstract Numerical investigation of natural convection within porous square enclosures has been performed for various thermal boundary conditions based on thermal aspect ratio on bottom and side walls. Penalty finite element analysis with bi-quadratic elements is used to solve the governing equations. The numerical solutions are studied in terms of streamlines, isotherms, heatlines, local and average Nusselt numbers for a wide range of parameters Da (10 −5 –10 1 ), Pr (0.015–1000) and Ra ( Ra = 10 3 –10 5 ). At low Darcy number ( Da = 10 −5 ), heatlines are perpendicular to the isotherms indicating conduction dominant heat transfer. As Da increases to 10 −3 , convection is initiated and the thermal mixing has been observed at the central regime for all A s. At low Prandtl number ( Pr = 0.015) with high Darcy number ( Da = 10 −2 and Da = 10 1 ), multiple circulations are observed in streamlines and heatlines and they suppressed for higher Prandtl number ( Pr = 1000). Isotherms are highly compressed along bottom wall at higher Prandtl numbers ( Pr = 0.7 and 1000) at A = 0.1 and 0.5. Temperature gradient is found to be high at the center of the bottom wall for A = 0.1 due to dense heatlines at that zone and that decreases as A increases from 0.1 to 0.9, irrespective of Pr , Da . Also, the temperature gradient is smaller at the top portion of side walls for A = 0.1 due to sparse heatlines along those zones and that is high for A = 0.9 due to dense heatlines. Distribution of heatlines illustrate that significant heat transport occurs from hot bottom wall to the top portion side walls at higher Darcy number ( Da = 10 1 ). It is found that Nu b attains maximum at X = 0.5 and minimum at corners for Da = 10 −5 , whereas that exhibits sinusoidal variation for Da = 10 −3 and Da = 10 1 irrespective of Pr and A . It is also found that Nu l follows wavy pattern at low Prandtl number ( Pr = 0.015) with higher Darcy number ( Da = 10 1 ) irrespective of A due to larger gradients of heatfunctions at several locations of left wall. The average Nusselt number show that the overall heat transfer rate is high at A = 0.1 compared to that of A = 0.5 and A = 0.9 irrespective of Da and Pr due to larger gradients of heatfunctions at A = 0.1.