3-D geochemical subsurface models, as constructed by spatial interpolation of drill-core assays, are valuable as-sets across multiple stages of the mineral industry's workflow. However, the accuracy of such models is limited by the spatial sparsity of the underlying drill-core, which samples only a small fraction of the subsurface. This lim-itation can be alleviated by integrating collocated 3-D models into the interpolation process, such as the 3-D rock property models produced by modern geophysical inversion procedures, provided that they are sufficiently re-solved and correlated with the interpolation target. While standard machine learning algorithms are capable of predicting the target property given these data, incorporating spatial autocorrelation and anisotropy in these models is often not possible. We propose a Gaussian process regression model for 3-D geochemical inter-polation, where custom kernels are introduced to integrate collocated 3-D rock property models while address-ing the trade-off between the spatial proximity of drill-cores and the similarities in their collocated rock properties, as well as the relative degree to which each supporting 3-D model contributes to interpolation. The proposed model was evaluated for 3-D modelling of Mg content in the Kevitsa Ni-Cu-PGE deposit based on drill-core analyses and four 3-D geophysical inversion models. Incorporating the inversion models improved the regression model's likelihood (relative to a purely spatial Gaussian process regression model) when evalu-ated at held-out test holes, but only for moderate spatial scales (100 m).