Incorporating endmember variability and spatial information into spectral unmixing analyses is important for producing accurate abundance estimates. However, most methods do not incorporate endmember variability with spatial regularization. This paper proposes a novel 2-step unmixing approach, which incorporates endmember variability and spatial information. In step 1, a probability distribution representing abundances is estimated by spectral unmixing within a multi-task Gaussian process framework (SUGP). In step 2, spatial information is incorporated into the probability distribution derived by SUGP through an a priori distribution derived from a Markov random field (MRF). The proposed method (SUGP-MRF) is different to the existing unmixing methods because it incorporates endmember variability and spatial information at separate steps in the analysis and automatically estimates parameters controlling the balance between the data fit and spatial smoothness. The performance of SUGP-MRF is compared with the existing unmixing methods using synthetic imagery with precisely known abundances and real hyperspectral imagery of rock samples. Results show that SUGP-MRF outperforms the existing methods and improves the accuracy of abundance estimates by incorporating spatial information.