In this work, a macroscopic nucleate pool boiling model was developed in OpenFOAM and validated against experimental data. This model takes into account the freeboard region, which is the region above the liquid free surface. The validated model was used to investigate the boiling physics in bulk liquid during pool boiling under surface heating (i.e. conventional heating). The simulation results showed that the turbulence dispersion force that acts on dispersed vapor bubbles can cause aeration in the bulk liquid, which can in turn influence the nature of the evaporation process. Besides, the results indicated that regardless of the amount of heat flux, the bulk liquid can barely be superheated during pool boiling under surface heating as vapor bubbles are generated exactly at the location of heat source (i.e. at heated wall). These bubbles immediately dissipate the superheat in the bulk liquid through evaporation cooling and thus prevents liquid superheating. It was also discovered that the evaporation at the liquid free surface can play a significant role in the overall phase change process, but it becomes less significant as the heat flux increases. Additionally, the roles of freeboard region in pool boiling simulation were assessed. The results showed that the inclusion of freeboard region in pool boiling simulation is necessary in order to achieve better mass conservation and provide a better prediction of the fluid flow.