In most modern coal mines, there are many coal quality parameters that are measured on samples taken from boreholes. These data are used to generate spatial models of the coal quality parameters, typically using inverse distance as an interpolation method. At the same time, downhole geophysical logging of numerous additional boreholes is used to measure various physical properties but no coal quality samples are taken. The work presented in this paper uses two of the most important coal quality variables-ash and volatile matter-and assesses the efficacy of using a number of geostatistical interpolation methods to improve the accuracy of the interpolated models, including the use of auxiliary variables from geophysical logs. A multivariate spatial statistical analysis of ash, volatile matter and several auxiliary variables is used to establish a co-regionalization model that relates all of the variables as manifestations of an underlying geological characteristic. A case study of a coal mine in Queensland, Australia, is used to compare the interpolation methods of inverse distance to ordinary kriging, universal kriging, co-kriging, regression kriging and kriging with an external drift. The relative merits of these six methods are compared using the mean error and the root mean square error as measures of bias and accuracy. The study demonstrates that there is significant opportunity to improve the estimations of coal quality when using kriging with an external drift. The results show that when using the depth of a sample as an external drift variable there is a significant improvement in the accuracy of estimation for volatile matter, and when using wireline density logs as the drift variable there is improvement in the estimation of the in situ ash. The economic benefit of these findings is that cheaper proxies for coal quality parameters can significantly increase data density and the quality of estimations.