A numerical investigation of heat conduction and laminar natural convection in ice-water systems containing porous metal foams, undertaken in the context of computationally convenient two-dimensional steady-state problems, is presented in this thesis. The overall goals of this work are to provide improvements to available cost-effective mathematical models of these phenomena, solve these models numerically, and investigate the influence of the porous metal foam on fluid flow and heat transfer in ice-water systems. The long-term goal (and the motivation for this work) is to contribute to the development of mathematical models and numerical solution methods for simulations of enhanced ice-water seasonal coldstorage systems.The proposed mathematical models are based on the local volume-averaging method. A Darcy-Brinkman-Forchheimer model is used for the momentum equations. For the heat transfer, volume-averaged equations governing two intrinsic phase-average temperature fields are used: one for the metal foam and the other for the water (solid or liquid). The following improvements to available two-temperature models are proposed: novel expressions for the interfacial heat transfer coefficient in both the conduction and convection regimes; and modified effective thermal conductivity models that provide consistency between predictions of one-temperature and two-temperature models in the limit of local thermal equilibrium.A well-established fixed-grid, co-located, finite volume method (FVM) is adapted for the numerical solution of the aforementioned mathematical models. All of the computer simulations are done with rectangular calculation domains, cooled and heated on the opposite side walls, and the adiabatic condition is imposed on the top and bottom walls.The FVM is first validated by the comparing the predicted results to experimental data for steady-state conduction and laminar natural convection in square enclosures containing pure liquid water and ice-water systems (no foam), with temperatures spanning the density inversion point of water. The problem involving natural convection in pure liquid water is solved using a variable-property model (VPM) and also a constant-property model (CPM), with the constant fluid properties evaluated at several reference (or average) temperatures,