A one-dimensional heat flux model has been developed in this study to analyze changes in heat propagation within debris cover on a glacier and the subsequent melt. The model can accept data collected in situ for calculations of real modeled melt, and also synthetic data based on simple wave forms and debris physical characteristics. The model is based on the Crank-Nicolson finite difference method of solving the one-dimensional heat equation. It has been written and implemented using MATLAB 2007b. A GUI has also been developed to allow for quick visualization and model adjustment. Numerous hypothetical runs have been performed to test the limitations of the model and to determine the accuracy and performance of the model. The model has been used to study what differences arise when using the melt calculation method of assuming an average daily linear thermal gradient versus using a physically based model. From this it has been shown that under stable atmospheric conditions and a relatively thick (~>0.5m) debris layer, the linear method does produce results similar to the full physical model. However, during times of unstable atmospheric conditions, errors do arise as the internal heat storage begins to fluctuate rapidly. If using the linear method on thinner (~<0.5m) debris layers the errors grow in response to the underestimation of internal temperatures. Using the model to analyze the changes in heat flux with varying debris thicknesses, a relationship between the two has been found, such that the linearity of the temperature profile is inversely related to the debris thickness. From this, a secondary relationship has also been found where the magnitude of the thermal gradient is also inversely related to the debris thickness. A comparison of calculated melt based on the assumption of a linear thermal gradient with that of the actual model shows a decreasing error with depth up to ~0.6m. It has been found that with an increasing debris thickness two things happen which decrease the error between the two methods. Firstly the magnitude of the modeled thermal gradient, at or near the surface, has been shown to decrease by almost a factor of 4 in this case, reducing the error between the two methods within this region of the debris. Secondly, the values for modeled thermal conductivity at the debris-ice interface approach zero, thereby reducing the error between the two calculated melts in the lower regions as the debris thickens.The model has also been used in an attempt to replicate results from a previous laboratory study performed by Reznichenko et al. (2010) where under stable conditions in a lab setting, melt was calculated under debris layers of varying thicknesses. The results of the model have been used as a verification of the accuracy of the model. Using data from the 2010 melt season from Longyearbreen, melt calculations have been made using both the modeled and the linear thermal gradient method producing values of 0.54m and 0.47m respectively. A correction has been made to the melt calculated using the linear thermal gradient based on an analysis of the differences between the two gradients and a correlation with daily average surface temperatures. By applying this correction equation the final R2 between the melt values was increased to 0.98. An Østrem curve has also been created using the model and varying the debris thickness by sub-centimeter increments. Based on the model’s performance tests and ability to produce similar melt results as laboratory results, the Østrem curve is believed to be more accurate in its general shape rather than its magnitude. The model has also been used to study the effects of measuring melt below varying debris thicknesses during shorter time periods throughout the melt season and how this subsequently changes the resulting Østrem curve. Differences on the order of ~0.5m have been found using 10 day periods during the mid-beginning melt season, mid-peak and mid-end of the melt season.