Phase change problems with multiple fronts in cylindrical systems, which arise in ice-on-pipe brine thermal energy storage systems, are investigated in this study. Two numerical methods, the boundary element method ( BEM ) and the thermal network method ( TNM ), are developed to solve the multiple front phase change problems. In the thermal network method, the distributed effects of sensible energy are approximated as lumped capacity at the boundaries and quasi-steady assumption is used. In the boundary element method, the full effects of sensible energy are precisely considered. The boundary element method is developed for multiple front phase change problems in one-dimensional radial systems. This method is applied as a module to a 2- dimensional axisymmetric problem and the model is used to predict the dynamic performance of the ice-on-pipe thermal storage systems with a parallel tube arrangement. The thermal network method is developed to solve the 2-dimensional axisymmetric problems with a moving external boundary and the model is used for the ice-on-pipe thermal systems with both counter flow and parallel flow arrangements. Performance predictions generated with the TNM are compared to experimental data from Oak Ridge National Laboratory (ORNL) for Calmac and Baltimore Air Coil (BAC) ice thermal storage systems which have counter flow arrangements. The predicted and measured results for the brine exit temperatures and load profiles are in good agreement throughout most of the charge and discharge cycles. Due to the lack of the experimental data for the counter flow arrangements, the results from BEM for the parallel arrangement problems are tested against the TNM and the negligible sensible heat approximation method (NSH) by comparing the outlet temperature, the latent state-of-charge (L-SOC) of the tube. It is found that BEM and TNM both give consistent and accurate results. / Ph. D.