Geostatistics is the most commonly used probabilistic approach for modeling earth systems, including quality parameters of various geoenergy resources. In geostatistics, estimates, either on a point or block support, are generated as a spatially-weighted average of surrounding samples. The optimal weights are determined through the stationary variogram model which accounts for the spatial structure of the samples. Recently, efficient modeling workflows using various machine learning algorithms (MLAs) have been expanded to the spatial context for modeling geological heterogeneity. The flexible use of MLAs as a spatial estimation tool stems mainly from the fact that unlike kriging, they do not require any variogram, nor do they depend strongly on a prior stationarity assumption (i.e., second order stationarity). This study evaluates the performance of two MLAs (ensemble super learner and elliptical radial basis neural network), ordinary kriging, and hybrid spatial modeling approaches using ordinary intrinsic collocated cokriging. The aforementioned modeling techniques are compared for estimating resources for four coal variables (wash yield, ash yield, calorific value and thickness) as an example. The results suggest that MLAs, when implemented alone, do not outperform ordinary kriging, but the estimation accuracy of the final model, measured by the root mean squared error tends to subtly improve (