Natural convection significantly affects the thermofluid transport processes for a material experiencing phase change. In the present work, steady-state solutions of the convection and associated momentum and energy transfer are computed by considering the solid, liquid, and mushy zones simultaneously. Finite-volume computational techniques, including the use of non-orthogonal curvilinear coordinates, second-order discretization, adequate number of nodal points, and adaptive grid solution method are utilized for wide ranges of Rayleigh, Prandtl, and Stefan numbers. It is found that the size and strength of the convection cell, as well as the location and shape of the phase boundaries are all strongly dependent on the combination of the above controlling parameters, as well as the temperature range governing the existence of the mushy zone. A systematic presentation of the effect of those factors on the transport characteristics is made to delineate the physical mechanisms responsible for the phase change process.